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method is investigated first. However, it is found that 
the results are heavily dependent on the location of the 
collocation points and characteristics of the equations. 
Therefore, the method is rejected as unreliable. Next, the 
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analysis of the pressure sensitive combustion instability, 
is considered. This method is found to work very well. It 
is found that the pressure wave forms exhibit a strong 
second harmonic distortion and a variety of behaviors are 
possible depending on the nature of the combustion process 
and the parametric values involved. Finally, a one- 
dimensional model provides further insight into the problem 
by allowing a comparison of Galerkin solutions with more 
exact finite-difference computations. 
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Chapter 1 


INTRODUCTION 

The steady operation of a liquid-propellant rocket 
engine is often disturbed by the occurrence of large 
pressure oscillations in the combustion chamber. These 
oscillations, which can lead to damage to or failure of the 
motor, are caused by amplification of initially small 
acoustic disturbances due to the energy released by unsteady 
burning of the propellant. This is usually referred to as 
combustion instability. It is generally accepted that 
unsteady burning can be correlated with the gas pressure 
(pressure sensitivity) and the magnitude of velocity of the 
fuel drops relative to the gas (velocity sensitivity). 

This dissertation is concerned with mathematical modeling 
of velocity-sensitive combustion instability in liquid- 
propellant rocket motors. 

A survey of literature dealing with mathematical 
modeling of pressure-sensitive combustion instability can 
be found in the dissertation of Powell (1). One of the 
first papers to examine nonlinear effects is that of Maslen 
and Moore (2) who considered only the fluid mechanical 
effects in a circular cylinder. The paper by Priem and 
Guentert (3) is one of the first to include the effect of 
velocity sensitivity. They discussed combustion instability 
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in a thin annulus. 

Since direct numerical solution of the governing 
equation^ is very time consuming and difficult to extract 
information from, many investigators such as Powell and 
Zinn (4), Lores and Zinn (5), and Peddiesion, Ventrice, and 
Purdy (6) have employed methods of weighted residuals such 
as Galerkin and collocation methods. 

Previous applications of the method of weighted 
residuals to combustion-instability problems have dealt with 
pressure-sensitive combustion. In this work this method will 
be extended to handle situations in which velocity sensitivity 
is important. Both the collocation and Galerkin methods will 
be considered. Stability boundaries and pressure wave forms 
will be computed numerically for transverse motion in a - 
cylindrical combustion chamber. In addition a finite- 
difference method will be used to solve the one-dimensional 
problem of transverse motion in a thin annular chamber. 

These results will be compared with those obtained by a 
method of weighted residuals solution of the same problem 
in order to assess the accuracy of the approximate solution. 

In this study, attempts to solve the velocity- 
sensitive combustion instability problem by various 
numerical methods are considered. The governing equations 
that describe the flow conditions inside liquid-propellant 
rocket motors will be derived in Chapter Two. An analytical 
solution is obtained in Chapter Three for nonlinear acoustic 
motion in a cylindrical chamber. The collocation method is 



applied to combustion instability in a cylindrical chamber 
in Chapter Four. In Chapter Five, the Galerkin method is 
applied to the same problem. In Chapter Six, attention is 
given to the analysis of combustion instability in a thin 
annulus. This allows the comparison of the Galerkin method 
and a finite-difference solution in a one-dimensioned 
context. Finally, a summary of this research is contained 
in Chapter Seven. 



Chapter 2 


COMBUSTOR EQUATIONS 


It is the objective of the discussion presented in 
this chapter to analyze the combustor conservation equations 
and reduce them to a tractable system. The simplification 
will be done in such a manner that will allow the resulting 
equations to retain both the mathematical and physical 
essence of the original problem. 

Under the assumptions that the usual balance 
principles of mechanics can be applied separately to each 
phase, the following equations are derived by applying the 
laws of conservation of mass, momentum, and energy to an 
arbitrary stationary control volume. 

Conservation of mass of the fluid requires that the 
time rate of change of mass in a volume v equals the rate 
at which mass enters v plus the rate at which mass is 
generated inside v and can be expressed as follows (for a 
fixed volume): 

* 

pdv - — f pu*nds + / wdv 
^ v s v 


-* * 
where n is an outward unit vector normal to s , p is the 

4 * 

density of fluid, u is the velocity of fluid, w is the rate 

* 

at which mass is generated per unit volume, and t is the time. 


ij 



Using the divergence theorm and the arbitrariness of 
the stationary control volume, yields 


* 

9 t* p 




( 2 . 1 ) 


Similarly, the balance law of mass for fuel phase is 

* * 

9^.*p^ + ?-(p L u L ) = - w (2.2) 

* : 

where p is the density of fuel and u is the velocity of 
L L 

fuel. 

The balance law of linear momentum states that the 
rate of change of linear momentum in the volume v equals 
the rate of entry of linear momentum across the surface s 
plus the sum of all external forces , acting on the volume v 
and can be expressed as 

** ** * * #* *+ 

9^#/ pudv = - / pu(u*n)ds + / ^dv + / wu T dv - / Pnds 
t* v s v v L s 

* 

where ? is the force per unit volume applied to the gas by 
x 

the fuel and P is the pressure of the gas . 

Using the same argument, the equation becomes 

****** # # 

9 t *(*u) + ^*(puu) = ? + wu L - VP. (2.3) 


Similarly, the balance of momentum for the fuel phase is 
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* * * * # „* 

5 t* Cp lV + ,,(p A S l ) = - 9 - ’S 


(2.4) 


Substituting (2.1) in (2.3) yields 


* ** 


p(3 t# u + u*^u) = P + w(u L 


# 

u) 


?P. 


(2.5) 


Similarly for the fuel phase 


* ** 


Pl ( 3 4-* u L + u L* Vu L } = 


* 

-> 

- F. 


( 2 . 6 ) 


The law of conservation of energy states that the 
rate of change of energy in a volume v equals the rate at 
which energy enters the volume v plus the rate at which 
energy is generated internally plus the rate at which work 
is done by external forces and gives 


% £ £ £ £ 

8.*/ p(e + ^u 2 )dv = - / p(e + ?iu 2 )(u*n)ds - / (Pn)*uds 
t * V s s 

# # * 

+ /^‘U^dv + -f v w(e L + %u£)dv + /^Qdv 

£ £ ^ ^ 

where e + ^u 2 is the energy of gas per unit volume, e + ^su 2 

Li Lt 

is the energy of fuel per unit volume, and Q is the heat 
transfer rate from the fuel to the gas. 


Following the same argument , the equation becomes 
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# * „# 

3 *Cp ( e + ^u 2 ) ) + V • (p (e + igu 2 )u) = -V • (Pu) 


+ $.u L + w(e^ + %u£) + Q. 


( 2 . 6 ) 


Similarly, the balance of energy for the fuel phase becomes 


3 


t 


* 




+ ^u 2 ) ) 
Li 


“ it Jt 

f - ( *L (e L 


+ ^U 2 )u ) 
L L 


*- * - 

+ * „ # 


P«u t — w(e_ + ^£U 2 ) — Q. 

Li Li Li 


(2.7) 


Substituting (2.1), (2.2) into (2.6), (2.7) respectively 
yields 


P *(3 #(* + ^u 2 ) + u»v(e + %u 2 )) = - v*(Pu) 
t* 


£ $ 

+ ^.u T + w((e + ?su 2 ) - (e + Jsu 2 )) + Q 
L L L 


( 2 , 8 ) 


# # 

P L C V (e L + %U l’ V* ( *L 


* * 


+ ^U 2 )) = 


- P-t - 5 . 

L 


(2.9) 


It is more convenient to work with the thermal energy 


equation for each phase. For gas phase this is done by 
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* 

dotting u into (2.3) to obtain the mechanical energy 
equation and subtracting this from (2.8) and using (2.1) to 
get 

* * * * * 

*De/Dt*= PD /Dt * ( log* ) + ?*(u T - u) + Q + w((e - e) 

Li L 

*2 *2 * * * * * 

+ (u L - u )/2 + u-(u - u L ) - P/p) • (2.10) 

In a similar way, it can be shown that 

" *: D T/D*t "(2 . Il)' 

L L L 

* * * # 

where D /Dt* = 9 *+ u*^, D T /Dt* = 3 + u r are the 

t* L t« L 

comoving derivatives for each phase. 

For simplicity, gas phase viscosity and heat 
conduction were neglected in the previous analysis. Since 
these produce dissipation, equations derived in this way 
should give a conservative estimate of the stability of the 
system. 

It appears (see, for instance, Powell (1)} that the 
primary effect of the burning fuel is that of interphase 
mass transfer. Thus, the balance laws for mass, momentum, 
and thermal energy for each phase will be further 
simplified by neglecting all interphase transfer terms 
other than those appearing in the mass-balance equations. 
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This leads to the system of equations : 


* * *11 * 
Dp/Dt + pv-u = w 


* * 

D* L /Dt* + P L ^ # ^ L 


* 

w 


* ** 

*Du/Dt * = - $P 

* * * 

p L^L^L^ Dt = 0 

*C V DT/Dt * = P(D /Dt *) (log* ) 

* * * 

p C DT /Dt = 0 (2.12) 

L 'L L 

* * * * 

where the constitutive equations e = C v T and e^ = 

. * 

have been used (with T denoting temperature and C y specific 
heat). The equation of state for an ideal gas is 

* * * 

P = pRT (2.13) 

where R is the gas constant. 

The governing equations will now be nondimensionalized 
with respect to a steady reference state, which will be 
denoted by the subscript "r" . All lengths will be referred 
to some characteristic length L r , such as the chamber radius. 
The characteristic velocity is the speed of sound at the 



reference state, and the characteristic time is the wave 
travel time L r /C^, The dimensionless quantities are 
defined below: 


* 


^ =v /L r 

* - Lpt/Cj. 

* 


II 

o 

S L“ C r S L 

# 

* 

p 3 p^p 

P L = p r p L 

# 


•n 

ii 

*-d 

w = p r C r w/L r 

* 2 

* 

T = C r T/( Y R) 

t L= °r T L / <f R > 

2 

2 

C r = YP r/ p r 

P p-pC-p/y • 


The dimensionless conservation equations become 
Continuity 

Dp/Dt + pV*u = w 

DpjVDt + p^'U^ = - w (2. 

Momentum 


Du/Dt = - ^P/ Y 
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Du L /Dt = 0 


(2.15) 


Energy 


DT/Dt = (y - 1)P (D /Dt ) ( logp ) 

DT L /Dt =0 (2.16) 


Equation of state 


P = pT • 


(2.17) 


By solving (2.l6a) one obtains 


T = pl-i . (2.18) 

Equation (2.17) then becomes 

P = p y • (2.19) 

Substituting (2.19) to (2.15a) produces 

Du/Dt = -V( p y_ 1 )/(y- 1) • (2.20) 

Taking the curl on both sides of equation (2.20) gives 

^ X Du/Dt = 0 • (2.21) 

It can be shown using vector identity that (2.21) 
leads to the following equation which describes the 
generation of the vorticity ft = ^ X u in the flow. 


Dft/Dt = fi-(v u) - fi( V-u) . 


( 2 . 22 ) 
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This equation describes the variation of vorticity 
for a given fluid particle. Suppose that at t = 0, this 
particle has zero vorticity and at the point where it is 
located the velocity gradients are bounded. Then by (2.22) 
the rate of change of vorticity of the particle vanishes. 

It follows from this that the vorticity will be zero at the 
next instant. By induction, it is seen that as long as the 
velocity gradients are bounded at each point occupied by the 
particle, its vorticity will remain zero for all time. If 
all fluid particles in the system have zero vorticity at 
t = 0, the vorticity will vanish at all points in the flow 
field and for all t >_0, therefore, as long as the velocity 
gradients are bounded through the flow field ( as they would 
not be, for example, at a shock wave ). Assuming this 
condition to be satisfied implies irrotational flow (fi = 0) 
and allows the definition of a velocity potential <f> defined 
such that 

u = V <}> . (2.23) 

Substituting (2.23) to (2.20) gives 

9 £<}> + +pl -1 / (y-1) = F(t) . 

Since the velocity is the space derivative of <f>, it 
is permissible to add to <f> any arbitrary function of time. 
This is equivalent to a statement that F(t) is arbitrary. 

For future convenience F(t) was selected to be l/(y-l). 

This choice yields 
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p Y ” i = l - CY-i)Ca t <c + %$♦•$♦) • (2.24) 

Substituting (2.24) into (2.14a), (2.18) and (2.19) gives 

3 d) - V 2 <|>(1 - (y- 1)(3 d> + *s7+*7*)) + V<t. • V (2 3 d> 

, u u u u 

+ V ^ • V 4> ) + (1 - (y-l)(3 t <t> + % V <J) * V <J> ) ) (y - 2)/(y-1) wb o 

T = 1 - (y-1)(3^4 i + ^<J»*^(|>) 

P = (1 - (Y-l)(» t * + ^<J>*^)) y/(y_1) • (2.25) 

This system of equations contains nonlinear gas. 
dynamics terms of all orders. 

Due to their highly nonlinear and mathematically 
complicated nature, the system of equations obtained above 
cannot be solved exactly. In order to obtain simpler, but 
approximate, equations which can be more easily analyzed, all 
nonlinear terms of order higher than two are neglected. 

This has the effect of retaining the most important 
nonlinear effect, while greatly simplifying the form of the 
equations. The system of equations then becomes 

p = 1 - (3 t <t> + + %(2 -y ) ( ) 2 

T = 1 - (y— 1)(3 d> + 

0 
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P - 1 -y( d t <f> + 4> * v <t> ) + hy(%^) 2 

3 tt 4>-V 2 4> (1 - (Y-l)3 t + ) + 2v^-vO t <J.) 

+ (1 - (y-2)3 a)w = 0 . (2.26) 

X- 

For a cylindrical combustor, it is desired to 
investigate the stability of the steady-state solution of 
(2.26d). To do this the steady-state solution must be 
obtained. It is assumed that the steady -state burning rate 
depends on z only (w = w(z)) and that the flow is axial 
( 4> = iCz)). Thus (2.26) simplifies to 

du/dz = w (u = d£/dz) 

T = 1 - hCy-l)u 2 
P = 1 - Jgyu 2 

£ = 1 - ku 2 . (2.27) 

It can be seen from (2.27) that the deviations of the 
steady-state solution from a uniform state are 0(u 2 ). To 
allow a discussion of orders of magnitude of various terms 
define the ordering parmeter e such that 


e = max (u) . 


( 2 . 28 ) 
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This Implies that 

u = 0(e) , w = 0(e), i = 0(e). (2.29) 

The stability analysis will now be carried out by assuming 
that 

= ±C z) + e<fi ’ ( x, y, z,t ) 

w = w(z) + e 2 w’ (x,y, z,t ) ( 2 . 30 ) 

where <f>' =0(l),w'=0(l). (2.31) 

It is being assumed that the unsteady perturbation of the 
velocity potential from the steady state is of the same 
order of magnitude as the deviation of the steady state 
from a uniform state and that the unsteady burning rate 
perturbation is of the same order as the unsteady nonlinear 
gas-dynamic terms. The first assumption Is necessary to 
be consistent with the quadratic approximation inherent in 
(2;. 26) while the second eliminates nonlinear terms involving 
products of w' and These assumptions are made for 

simplicity and are not meant to imply that other orders of 
magnitude for the various terms are impossible or even 
unlikely. The objective of this work is to pick one case 
which is mathematically tractable and subject it to extensive 
analysis. While it is felt that most of the features of the 
stability problem will be reflected by this case, it is 
recognized that many other possibilities remain to be 
explored . 
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Substituting (2.30) and 

P=£ + ep ' , P = P + eP» , T=T+eT' (2.32) 

{ p * = 0(1), P* = 0(1), T' = 0(1)} into (2 . 26 ), retaining all 
terms of 0(e 3 ) or lower, and dividing all resulting equations 
by e leads to 

a tt r - V 2 $» + 2ud zt 4' + w a t <J) ' 

+ e { ( y-1 ) V 2 <J> 1 9 t <f>' + 2 V <f> ’ • V ( 9^. <}> ' ) + w' } = 0 
p’ = - + u8 z <(,’ + % e( V4> » - V4, » - (2-y) ( 9 t <(. ' ) 2 ) } 

p« = _ Y { 9 d> 1 + u 9 + he(V<p' • V4> * - (9 d’) 2 )} 

t — 2 t 

T* = -Cy-DO^' + U3^' + 3se7*' (2.33) 

The equations derived from this perturbation analysis 
are expected to be valid as long as the amplitudes of the 
perturbation quantities are finite and smaller than unity. 

For simplicity the primes appearing in the above 
equations will be dropped and it will be understood that 
unprimed quantities are associated with perturbation from 
the steady state. 



Chapter 3 


ANALYTICAL SOLUTION OP WAVE EQUATION 

In this chapter a second-order solution for the 
velocity potential associated with transverse acoustic wave 
motion of an ideal gas in a circular cylinder is obtained. 
The solution is found in the form of a Fourier-Bessel 
series. This provides a standard to compare the proposed 
methods which will be discussed in the following chapters. 

The governing equation is obtained from (2.33a) by 
assumption that the burning rate function is absent. Then 
the equation reduces to the second-order form of the one 
discussed by Maslen and Moore (2) for pure gas motion in a 
cylinder. 


— V^<J> + eC2v<|> • 7 9^. 4> + (y— 1) V^<t> 9 ^ ) = 0 (3*1) 

For a cylinder of radius L r , a set of cylindrical 
polar coordinates (r,e,z) with the z axis being coincident 
with the cylinder’s axis of symmetry is used. In seeking 
a second-order transverse wave solution, one assumes 

<J> = <J> (r,e,t) + e<f> (r,e,t) + ... . (3.2) 

1 2 

Substituting C3.2) into (3*1) and (2.33) leads to 
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a tt«, - ,2 *, = 

1 

3 t t* - V 2 * = - (a f (u 2 ) + (y-1)72* 3 * ) (3.3) 

^ 2 2 c 1 1 r 1 

P = 1 + eP + e 2 P (3.4) 

1 2 

where P 2 « - y3 * , P = -y<M + ^(u 2 (3 <{. ) 2 )) (3-5) 

tl 2 * 2 1 t 1 

72 ^i = 3 rr+i + 3 r*i /r + 9 ee <J>j/ r2 > 1 = 1 » 2 

u 2 = + t 3 ©^/ 1 ') 2 . (3.6) 

Three different cases are discussed below. For the 
initial conditions 

*(r,e,0) = J (S r )cos0 , 3.*(r,0,O) = 0 (3-7) 

ill t 

the solution of C3.3a) is 

<f> = J (S r)cos0cos(S t) 

l i n li 

Substituting this equation to (3* 3b) produces 

3 t1 .<j> — V 2 4 > =(b + z z J m (S mn r)cos(m0))sin(2S t) (3.8) 

2 2 oo m= q n= i m mn li 

m^i 

where J is used to denote an n’th order Bessel function of 
n 

the first kind, S mn is used to denote the m’th root of 
the equation JH = 0, and a prime indicates differentiation 
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with respect to the argument of the Bessel finction. The 

coefficients b are integrals of Bessel functions which 
mn 

must be computed numerically. The details of this 
calcuation are discussed in appendix A. 

Solving C 3 - 8 ) one obtains 

♦ = (b /(4S 2 ))(2S t - sin(2S t)) 

2 oo 11 11 11 

+ z ? (b / (S 2 - 4S 2 ))J (S r)cos(me)Csin(2S t) 

q 2 mn mn 2 2 m mn i i 

m / 1 

- ( 2 S 1 / S mn )sln(s mn t)) - (3 ' 9) 

The other two solutions to be discussed subsequently are 
found in the same way. For the initial conditions 

/ i 

<t>(r,6,0) = 0, ^<f>(r,9,0) = S J (S r)cos0 (3-10) 
t T 11 l 11 

it is found that 

<p = J (S r)cos0sin(S t) (3.11) 

l l 11 11 

and that <p 2 is given by the negative of (3*9)- For the 
initial conditions 

<p(r,6,0) = J (^ 12 r )cos0 j 

(r , 0 , 0 ) = S J (S r)sin0 (3*12) 

x 11 l 11 


the solution is 
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$]_ = J-^Cs^^jcosCs^-^t -e) 

and *2 = n li {2b 2n /(S 2n “^ll 5 }J 2 (S 2n r ) (sin(29 ) {cos (25^ ) 

- cos(S 2n t)} + cos(20){sin(2S l;L t) - (2S 11 /S 2n )Sin(S 2n t ) } ). 

It can be seen that for these initial conditions <j>j 
represents a spinning wave but <t> x + e<f> 2 does not. Some 
typical results for the pressure functions P-^ and P 2 were 
computed using the solutions for initial conditions (3-7) 
(labeled standing wave) and initial conditions (3-12) 

(labeled traveling wave) for r = 1, 0 = 0, and y = 1.2. 

These data are presented graphically in Figure 1. For 
these computations the series expansions were terminated at 
n = 5- There is virtually no difference between the results 
for n = 4 and n = 5, and even the results for n = 1 provide 
a reasonably accurate solution. 

The results discussed above were used to check the 
numerical calculations to be discussed in the next two chap- 
ters. These results generalize those of Maslen and Moore (2) 
to initial conditions which do not lead to periodic solutions. 



Chapter 4 


COLLOCATION 

In previous investigations of velocity-sensitive 
combustion instability {see, for instance, (6) and the 
references therein) the most widely used burning-rate 
equation is that associated with the vaporization limited 
combustion model discussed by Priem and Guentert (3)* In 
the notation of the present thesis the second-order version 
of this equation can be expressed as 

w = (w/ e 2 ) C (1 - %e3 t <J>)((l - u z /u L ) 2 +(u t /u L ) 2 ) iG -l) (4.1) 

where u is the component of perturbed gas velocity in the 
z 

axial direction and u is the perturbed component . normal to 
the axis of symmetry, and u L is the magnitude of the droplet 
velocity vector which assumed to be axial. Substituting 
(4.1) into (2.33a) yields 

+ + 2ud^ip - v 2 <f> ( 1 - e(y-l) ) 

+ 2£^(j)*^a ' <ji + (w/e)((l - 3 <J>)((1 - 3 <{>/u T ) 2 

t — t z u 

+ C v<f> • v<j> - (a 4 , ) 2 )/ u 2 )^_ l) = o. (4.2) 

Z 
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No exact solutions are to be expected to (4.2). Fully 
numerical solutions, on the other hand, are expected to be 
quite time consuming for multidimensional problems. For 
these reasons it was decided to attempt approximate 
partially analytical solutions using the method of weighted 
residuals {see for instance, (6)}. For problems of pressure 
sensitive combustion instability Powell (1) successfully 
employed the Galerkin method. This method is limited in 
applicability, however, to equations containing nonlinear- 
ities involving only polynomial functions of the dependent 
variable and its derivatives. It is clear that (4.2) is not 
of this form. The orthogonal collocation method was, 
therefore, chosen because of its simplicity and adaptability 
to equations containing nonlinearities of any algebraic 
form. The application of this method to the problem of 
transverse wave motion in a cylindrical chamber will be 
discussed below. 

Consider the transverse wave motion in a cylindrical 
combustor. It is convenient to describe this problem in 
terms of cylindrical polar coordinates (r,e,z). Then 

<t> = <p(r,e ,t) . (4.3) 

It will also be assumed that w and u L are constants. Then 
(3*2) becomes 

8 <p - + w 3 <|> — ( 3 <f> + 3 <p/r + 3 4 >/r^)Cl — e (y— 1) <)> ) 

tt — t rr p 0 0 t 
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+ 2e(a + 3„<j> 3 q . <f>/r 2 ) + (w/e ) ( (1 - ht d±. $ ) ((1 

r rt 0 0t — z 

+ C(8 r <|>) 2 + Ca e <J»/r ) 2 )/u 2 )* - 1) = 0 . (4.4) 


Now, assume an approximate solution of the form 


P Q 

<t> = I I J (S mvi r)cos(rn0)f (t). 
m=on=i ro nin mn 


(4.5) 


Equation (4.5) is a superposition of the normal modes 
associated with the corresponding linear acoustic problem. 

A solution of this form allows one to investigate the 
influence of nonlinearities on the behavior of modal 
amplitudes which would exhibit simple harmonic motion in 
the linear case. Only standing modes can be described by 
(4.5). To represent spinning modes another series involving 
sin (me) must be added. This is omitted for simplicity in 
this discussion. It is convenient to rewrite (4.5) as 


* 


N 

l J 


(S r )cos (m, 6 )f . (t) 

m j n j 3 3 


(4.6) 


where N = (P + 1)Q . 

The equation is now evaluated at N points (the collocation 
points) to yield 


*1 ' 


(4.7) 



where 


G . . = J (S_ „ r., )cos (m., e., ) . 
I J rrij m j ru l J 1 


(4.8) 


Next, (4.7) is inverted to obtain 


N -1 

f , = £ C. U. 

i j=i ij J 


(4.9) 


Then the various derivatives of (4.6) can be expressed in 
terms of value of <j>^ as 


N 


N 


(a r* J i i = i °ij ^ j 3 ^rr^i ili C ij*j 


N 


N 


C 3 e <|) ) ± = ^ 3 © e * ) i = 1 ii C ij <f> j i=1 »--- N CM- 10) 


where 


°L- = ^ilc^* = (8 r 1 r 1 C ik )C iij 1 


c ij (a e c ik )c kj J c ij = (a e i 0 1 c ik )c kj' 


(4.11) 


Substituting into (4.4), a set of N ordinary second-order 
differential equations having the following form is 
generated. 


N 


N N 


<f>. + w$, - £ C. , 4 . .(l-e(y-l)<j), ) + 2e £ £ C. 

i — i j = i ij J J j = lk=l 1 J K J k 


N N 


+ (w/e)( (l-he<p , ) (1 + Z Z C * * /u*) *-l)=0 (4.12) 

- l ,1 = 1 k=i ijk j k h 


where 


C V , r rr c r , C 00 /r2 

c ij C ij °ij /r i L ij /r i 


ij' 
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c ljk - c lJ c ik + °l)°ltc/ r l • (4 - 1 3 ) 

Solutions were obtained by a fourth-order Runge-Kutta 
method to the set of equations (4.12). It was found that 
the results were very sensitive to the number and locations 
of collocation points, especially to the radial distribution 
of points. It was found that even in cases when w was equated 
to zero (no combustion) it was possible to find many choices 
of collocation points which lead to the computed modal 
amplitudes increasing without bound. 

In order to illustrate the difficulties involved in 
a simple case, consider a one-dimensional problem of 
longitudinal wave motion with y = 1, u = 0 and w = 0. In 
this case (4.2) reduces to 

9 <t>-9 + 2e(3 4>a «j») = .0 . (4.14) 

tt zz z zt 

> 

First consider the boundary condition 

9 <|>C0,t) = 0 , <Kir,t) = 0 . (4.15) 

Zi 

A one-term solution satisfying (4.15) is 

<j> = f (t)cos( J iz) . (4.16) 

Applying the collocation method with a collocation point z 0 
produces an equation for <J> 0 = 4>(z 0 ,t) of the form 
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$+*$$+ ^(tan^Cz /2 )<p <t> ) = 0. (4.17) 

oo ooo 

The Galerkln method (to be discussed in detail in the next 
chapter) when applied to the same problem yields an equation 
for f of the form 


f + + 4 ff/(3ir) = 0. 

As another example consider the boundary conditions 

Sg^COjt) = 0 ) Bgfj^iTjt) = 0. 

A solution satisfying (4.19) is 

<t> = f(t)cos(z) 

for which the collocation method leads to 

.. 

<f> + d> - 2etan 2 z 4> 4> =0 (4.21) 

o o ooo 

and the Galerkin method leads to 

f + f = 0. (4.22) 

From the above two cases, one can observe that the 
last term in each of the equations obtained by collocation 
can take on any value between 0 and 00 depending on the 
location z 0 of the collocation point. No such difficulty 
is encountered when using the Galerkin method. Even if more 
points were used, the same behavior was obtained. Therefore 


(4.18) 


(4.19) 


(4. .20) 
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there appears to be no way to rationally decide which sets 
of collocation points can produce the correct result. For 
this reason, after the expenditure of much effort, the 
collocation method was abandoned and the Galerkin method 
was adopted. This will be discussed in the next chapter. 



Chapter 5 


THE GALERKIN METHOD 

The Galerkin method is a special application of the 
method of weighted residuals (usually referred to as MWR). 

It has been extensively used in the solution of various 
stability and aeroelasticity problems (see for instance 
Powell (1) } and proved Itself as a useful tool for the 
solution of both linear and nonlinear problems. Although 
it is an approximate mathematical technique, it has 
nevertheless produced results which were In excellent 
agreement with available exact solutions. These approximate 
solutions are usually simpler in form than the exact 
solutions obtained by numerical integration, and their 
quantitative evaluation requires considerably less 
computation time. However, this method is known to be 
reliable and applied conveniently only to equations 
involving nonlinearities of a polynomial type. Therefore, 
if this method is used in conjunction with the varporizat ion- 
limited burning rate function of (4.1) the problem becomes 
intractable. Thus, the following simpler purely phenomenolog- 
ical burning-rate function was employed for the case of 
instantaneous combustion response. 

w = wnu 2 (5.1) 
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where n is a constant which will subsequently be referred to 
as the interaction coefficient. It can be determined either 
by comparison of predictions of the theory directly with 
experiment or by comparison with burning rate laws meant to 
apply to special types of combustion processes. As an 
example of the latter method consider the vaporization- 
limited burning rate law (4.1). As it stands (4.1) exhibits 
both pressure and velocity sensitivity. A purely velocity 
sensitive law can be obtained by assuming e « 1 to get 

w v = ( w/ e 2 ) ( ( 1 + Cu/u l ) 2 ) H - i) (5.2) 

where the v is used to denote the vaporization model. If one 
equates the slopes of (5.1) and (5.2) at u 2 = 0 and plots 
both functions, the graph would look as follows. 
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It can be seen that (5.1) would overestimate the burning 
rate. Thus the upper bound for n can be computed. Prom 
(5.1) and (5.2) 

d w = wn, d ,W = (w/(4 e 2 u T 2 ) )/(l+(u/u T ) 2 ) 

U 2 — U 2 V — Li L 

Thus d ,w(0) = wn, d ,w (0) = w/(4e 2 u£). 

— U 2 V — Li 


Equating these results, one gets 


n „ = l/(4e 2 u 2 ). 

max L 


(5.3) 


The phenomonological law (5.1) can be related to other 
special burning-rate laws in a similar manner. 

If it is desired to consider history-dependent 
combustion processes , this can be done through the general 
formula 

t 

w = n/ G ( t — £ ) d (u 2 )d £ (5.4) 

0 5 

where G is memory function and £ is a dummy variable. If 
G(t) = H(t) (H being the unit step function) all increments 
of change in u 2 occurring in the past are counted equally 
and (5.1) is recovered. If G(t) = H(t) - H(t-x) all 
increments of change in u 2 are counted equally up to t units 
of time in the past while those previous to that time are 
not counted at all. Substituting this expression into (5.4) 



31 


and integrating one obtains 

w = wn(u 2 - u 2 ) (5*5) 

— T 

where u T =u(t-x ) and t is called the time delay. Equation 
(5*5) will be used to account for the history of the burning 
process in a rough way. More sophisticated treatments of 
this phenomenon could be obtained by employing a more 
realistic function for G(t). Substituting (5.5) and (5-1) 
into (2.33a) one obtains 

3 tt <|> - V 2 <}> + 2u3 zt <|> + W3 t <|> + e( ( y-1 ) V 2 <|> 9 t $ + 2V<t>*V3 t <|> 

+ wn(y<(>*7<f> - jV<J> ‘V<p )) = 0 ( 5 . 6 ) 

— XT 

where j = 0 for instantaneous combustion and j = 1 for 
history-dependent combustion. 

The most general solution of (5.6) (subject to hard- 
wall boundary conditions for the unsteady variables) can be 
written in the form of the following Fourier-Bessel series. 

+ = f 00^ + ) J 0^ S 0n r ^ 

+ £ e (f (t)cos(me)+g (t)sin(me))J (S r) (5-7) 

\ mn mn m mn 

Substituting (5-7) into (5*6) and applying the usual 
Galerkin orthogonalization procedure leads to an infinite 
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set of coupled ordinary differential equations governing the 

functions f and g . No solution can be obtained unless 

mn mn 

the series (5.7) is truncated so as to produce a finite set 
of equations. If the terms neglected are actually small, an 
accurate solution will result. 

In this thesis attention will be focused on initial 
disturbances having the form of the first tangential mode. 
The simplest finite series contained in (5.7) capable of 
modeling the effect of quadradic nonlinearities in (5.6) is 


♦ = fo (t)J o CS oi r) + f ) J 1 ( s 11 r ) cos9 + f 2 (t )J 2 (.S 2i r )cos26 
+ g 2 (t )J 1 (S 11 r )sin0 + g 2 (t )J 2 (S 21 r)sin20 • (5.8) 

The Galerkin method then produces the following ordinary 
differential equations governing these five variables. 

V a'o + s oi f o + eCc iVo + c 2 (f A + «A> 

+ c 3 C f 2 f 2 + g 2 g 2 ) + wntc^f* + c 12 (q + e{) + c 13 (f| + g|) 

- J<0 ll f OT + C 12 tf lx + Six' + C 13 <f Ix + Six”" ■ 0 


? 1 + * f l + S il f l + ‘“Wl + Vl f 0 + °6 (s i s 2 + f lV 


+ V g 2 S l + W + S nCC l4 f 0 f l + C 15 (f l f 2 + S 1 S 2 ) 
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- j( c i4 f 0 T f lT + C 15 (f lT f 2 T + s 1t s 2 t )))) " 0 


f 2 + - f 2 + S 21 f 2 + e ^ C 8 f 0 f 2 + C 9 f 2 f 0 + C 10^ S 1 S 1 f l f l^ 
+ HP( c 1 6^ f l“ S l^ +C 17 f 0 f 2 _ J '^ C l6^ f lT -S lT^ + C 17 f 0r f 2r^^ = ° 


S 1 + - s l + S ll s l + e ^ C 4 f O s l + C 5 S l f 0 + C 6^ f l s 2 g l f 2^ 


+ C 7 (g 2 f 1 f 2 S l^ + - n ^ C l4 f 0 S l + C 15^ f l S 2 S l f 2^ 


0^ C l4 f O T g i T + C 15^ f lT S 2 T g lx f 2x^^ ° 


g 2 + — g 2 + S 2 2 l g 2 + £CC 8 f 0 g 2 + C 9 S 2 f 0 * C 10 (g l f l 


+ f g ) + wn(C f g + 2C ,f g 

1 6 1 ~ 17 0 e 2 16 l e l 


- J(C 17V S 2 T 


2C l6 f l T s lx )))) ■ ° 


(5.9) 


The coefficients (1=1,2, .. .17) are integrals of 
Bessel functions which must be computed numerically. The 
details of this calculation are discussed in appendix A. 

The stability boundaries in the (n,e) plane 
were determined first. For a given value of w, a value of 
e was selected and solutions were obtained for various values 
of n. In each case, if the modal amplitudes exhibited 
growth with time, the value of n was decreased for the next 
run while if decay of modal amplitudes was observed, the 
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value of n was increased accordingly. A systematic 
iteration process was devised so that the computer could 
carry out these calculations without analyst intervention. 

In this way, one point on the stability boundary was 
established. Then a new value of e was chosen and the 
iteration process was repeated to determine another point 
on the stability boundary. This procedure (which consumed 
a considerable amount of computer time) was continued until 
enough points had been found to establish the shape of the 
stability boundary. It should be pointed out that the 
amplitude of the IT mode -always initially decreases due to 
the fact that purely velocity-sensitive combustion 
instability is always linearly stable. Thus it is necessary 
to continue the calculation for a considerable period of 
time to determine whether a given set of conditions 
corresponds to nonlinear stability or instability. 

Figures 2 and 3 show some typical stability 
boundaries for instantaneous combustion using the initial 
conditions 

f 1 (0) = 1, f 0 (0 ) = f 2 (0 ) = g 1 (0) = g 2 (0) = 0 
f 0 (0) = fiCO) = f 2 (0) = g^O) = g 2 (0) = 0 (5.10) 

and f x (0) = 1, f 0 (0 ) = f ^ (0 ) = g^O) = g 2 (0) = 0 

g x (0) = 1, f Q (0 ) = f (0) = f 2 (0 ) = g 2 (0) = 0, (5.11) 
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respectively. The region below and to the left of the 
stability boundary is associated with stable conditions, 
while the region above and to the right is associated with 
unstable conditions. 

In the case of linear acoustics, the first set of 
initial conditions would lead to a standing wave and the 
second set of initial conditions corresponds to a traveling 
wave. In both cases, it can be seen that the stability 
boundaries have roughly the forms of rectangular hyperbolas . 
As expected, increasing the steady-state burning rate reduces 
the value of n required to produce, instability. This effect 
is somewhat more pronounced for the traveling waves than for 
the standing waves. There does not appear to be a distinct 
pattern to the results. Thus for w = 0.2 standing waves are 
more stable than traveling waves while for w = 0.1 traveling 
waves are more stable than standing waves. 

Figures 4 and 5 show stability boundaries associated 
with initial conditions (5*10) and (5.11) , respectively, for 
history-dependent combustion. It is apparent that the 
influence of the time delay parameter on the results for 
standing waves is much greater than on the results for 
traveling waves. Again, no clear pattern emerges from the 
results. Comparing Figures 2 and 4, and 3 and 5 shows that 
instantaneous combustion can be either more or less stable 
than history dependent combustion depending on the value of 
the time-delay parameter. Comparing Figures 4 and 5, it 
can be seen that standing waves can be either more or less 
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stable than standing waves for history dependent combustion. 
Figures 4 and 5 also illustrate the fact that increasing the 
amount of time delay can either increase or decrease the 
stability of the system. Thus , for both standing and 
traveling waves, x = it corresponds to a more stable situation 
than either t = 0.5n or t = 1.5ir . The lack of patterns 
exhibited by these results emphasizes the need for numerical 
methods of the type developed during the present research. 

It appears that the only way to find out what will happen in 
a given situation is to solve the equations for that 
particular case. This matter will be discussed further 
subsequently . 

The pressure perturbation can be calculated from the 
velocity potential perturbation by substituting equation 
(5.8) into (2.33c), expanding the result in a Fourier-Bessel 
series and retaining only the terms corresponding to the IT, 
2T and 1R modes to obtain 


P ■ - > (( Vo + £(d 02 f 0 


d (f 2 + g 2 ) 

03 1 S l 


+ d 04 (f ! + S 2 } 


d 05^0 


+ d 06^ f l + + d 07^ f 2 


&2 ) ) ) Jq ^01 r ^ 


^ d ll f l + e ^ d 12 f 0 f l + d 13^ f l f 2 + s l s 2^ + d l4 f 0 f l 


+ d 15 Cf 1 f 2 + g 1 g 2 )))cos0 
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+ ( d nS;i + e ^ d i2 f 0 S l + d 13^ f> l S l ” S l f 2^ + d l4 f 0 S l 
+ d 15 ( f 1 g 2 - f 2 g 1 ))) sin0 ) J 1 ( S 11 r ^ 

+ (( d 21 ^ 2 + e ^ d 22 f 0 f 2 + d 23^ f l “ S l^ + d 24 f 0 f 2 

+ d 25 ^j " g£)) )cos2e 

• • • 

^ ^ d 2]_g]_ ^( d 22 ^*Qg 2 ^ d 23^i®i ^24^0^2 

+ 2d 25 f 1 g 1 ))sin20)J 2 (S 21 r)). (5.12) 

The coefficients d,. , d, , d„ are calculated in the same 

Os Is* 2s 

fashion as those discussed previously and the values are 
listed in table 8 of the appendix. 

Some typical results for wall pressure waveforms are 
presented in figures 6 through 37* 

Figures 6 through 9 correspond to instantaneous 
combustion with e = .05, w = .1, n = 175, and initial 
conditions (5.10). This leads to a stable standing wave 
oscillation and the pressure can be observed to decrease 
gradually. However, by changing the interaction index to 
n = 220 (Figures 10-13) the pressure is seen to grow 
gradually. This is an unstable situation. In both cases 
the response is dominated by the IT mode but distorted to 
some extent by the presence of the 2T and 1R components. 

Figures 14 through 17 correspond to history-dependent 
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combustion with n = 3^5, while the other parameters are the 
same as before. This is a stable situation and the amplitude 
of the pressure slowly decays with time. 

Figures 18 through 21 are obtained by changing the 
interaction index to n = 352. This is an unstable situation 
as indicated by the growth of the pressure amplitude with 
time. In both of these cases the presence of the 2T and 1R 
components is much more noticable than it was in the first 
two situations. 

Figures 22 through 25 correspond to instantaneous 
combustion with n = 180, e = .05* w = .1, and Initial 
condition (5.11). This produces a traveling wave. The 
pressure is observed to be decreasing with time (a stable 
situation) while changing the interaction index n to 210 
(Figures 26 through 29) causes the pressure to increase 
(an unstable situation). Figures 30 through 37 show the 
significance of the time delay function. Figures 30 through 
33 are obtained by setting r = ir and n = 197. This is a 
stable case. The pressure is observed to decrease. Setting 
n = 199 (Figures 3^ through 37), on the other hand, produces 
an unstable situation when the pressure increases very 
rapidly. In both of these situations the response appears 
to be dominated by 2T contribution to the pressure. 

From these figures one can conclude that the pressure 
waveforms exhibit a strong second harmonic distortion and 
this distortion arises from the effect of the quadratic 
nonlinear terms. A variety of behaviors are possible 



39 . 


depending on the nature of the combustion process and the 
parametric values involved. 



Chapter 6 


ONE-DIMENSIONAL MODEL 

In this chapter an annular combustion chamber with a 
gap width much smaller than the inner radius is considered. 
This geometry will henceforth be referred to as that of a 
narrow annulus. While this geometry is not of much practical 
interest, it is quite useful for the purposes of analysis. 

This is because only one space coordinate is needed to 
describe the problem. This makes a direct finite-difference 
numerical solution of the original partial differential 
equation feasible and also simplifies the algebra required 
to carry out the Galerkin modal analysis. In what follows 
three questions will be investigated. First, the effect of 
changing the numerical values of certain coefficients 
appearing in the governing equations for the modal amplitudes 
will be discussed. Second, the Galerkin solution will be 
checked using a finite difference numerical solution of the 
complete equation. Third, numerical solutions of the 
complete wave equation using the vaporization-limited 
burning-rate law will be compared to similar solutions 
associated with the phenomenological burning-rate law 
employed in the previous chapter. 

The appropriate wave equation for transverse wave 
motion in a narrow annulus can be obtained from (5.6) by 
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setting r = 1 and 3^ = 0. To further simplify the results 
the parametic values y = 1 (isothermal process) and j = 0 
(instantenous burning response) will be employed. For this 
situation (5*6) simplifies to 

a tt * + w3 t * - 9 0e 4> + 2e3 0 <(.8 et <). + wne(3 0 *) 2 = 0. (6.1) 

To carry out the modal analysis, it is assumed that 
the potential function can be expressed as 

i 

<p = f 1 (t)cos0 + f 2 (t)cos20 + g^(t)sin0 +g 2 (t)sin20. (6.2) 

Then applying the usual Galerkin orthogonalization 
procedure leads to 

+ wf x + eA 1 (f 1 f 2 + f 2 f 1 + g x g 2 + g 2 g a ) 

+ ewnA^f^ + g^g ) = 0 
• • • * • 

f 2 + «|f 2 + wfg + eA-(gjg^ - f 1 f 1 ) + ewnA lt (g^ - f|) = 0 

®1 + fi l S l + - S 1 + eA i (g 2^l + f l S 2 “ S 1^2 ~ f 2 S l ) 

+ ewnA 2 (f 1 g 2 - g 1 f 2 ) = 0 
• • * • • 

g 2 + Q^g 2 + wg 2 ~ eA 3 (g 1 f 1 + f^) - 2ewnA lt f 1 g 1 =0 (6.3) 
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where 8^ = l, 8^ = 2, A^ = 2, k^ = 2, A^ = 1 and A^ = .5. 

The symbols 8^, 8^, A^, A 2 , A^ and A^ have been inserted to 
illustrate the effect of changing their numerical values. 

McDonald (7) in his analysis of combustion- 
instability in. an annulus found that for a given value of e 
the value of n required to produce instability was approx- 
imately twice as high for standing waves as for traveling 
waves. The results discussed in the previous chapter for 
a full cylinder indicated no such relationship. Equations 
( 6 . 3 ) were employed in an attempt to determine the factors 
which have a significant effect on the stability boundary. 
Several calculations were made and a representative sample 
of the data thus obtained is presented in Tables 1 and 2. 

It was determined by McDonald that the terms representing 
gas-dynamic nonlinearities, had a small qualitative effect 
on stability calculations. Thus the coefficients A^^ and A^ 
were held fixed during these calculations. 

The entries in the first two lines of each table 
were computed using the correct equation for an annulus. It 
can be seen that the standing wave is twice as stable as the 
traveling wave, in agreement with the results of McDonald. 
The third and fourth lines in each table were computed by 
changing A^ from .5 to .77 (This makes the ratio A^/A^ 
the same as the ratio of the corresponding terms in the 
governing equations for the full cylinder.). It can be 
seen that this lowers the stability limit in all cases but 
does not alter the fact that an initial disturbance in the 
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form of a standing wave is twice as stable as one in the 
form of a traveling wave. 

The fifth and sixth lines in Tables 1 and 2 are found 
by restoring to its original value .5 and changing 
from 2.00 to 1.66 (This makes the ratio e< ^ ua -'- to the 

ratio associated with the full cylinder,). For this 

situation it can be seen that the standing wave is approx- 
imately ten times ais Stable as the traveling wave. Further- 
more, the effect of the value of ^2 on the location of 
stability boundary is much greater for a standing-wave 
motion than for a traveling-wave motion. The last two lines 
are associated with implementing both of the changes dis- 
cussed above simultaneously. These results confirm that the 
change in lowers the stability limit in all cases but 
does not affect the relative stability of standing-wave and 
traveling-wave disturbances. 

The data presented above indicate that standing waves 
will be twice as stable as traveling waves only under very 
special circumstances 3 2). There is no reason to 

expect this to be a characteristic of other systems exhibit- 
ing combustion instability as, in fact, it is not for a full 
cylinder. 

In this section, a finite different procedure is 
employed to generate a set of second order differential 
equations. The central difference formulas 


Vi 



*i-l )/(2A0) 
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3 0 e 4> i = ($ 1+1 - 2^ +$ i _ 1 )/(A0) 2 (6.4) 

are used to produce the following set of governing equations 
= -wi}>^ + (tf>^ + 2 “ + ” e (y-1) )/( A0 ) 2 

- e( *i + l - *i-l )( *i + l - *i-l )/(2Ae)2 

-ew i 2 < i < N - 1 (6.5) 

where N is the number of points employed. A forward 
difference formula 


9 0 <J>1 = (~34> 1 + 4<j> 2 — <J) ^ ) / (2A6 ) 


( 6 . 6 ) 


and backward difference formula 


~ ^N-2 ~ ^N-l + 3<J»j y j)/(2A0 ) 


(6.7) 


are used for the boundary conditions 


3 ft 4> CO ) = 0, 


3 o (u ) - 0. 


Hence, along the boundaries, (6.5) becomes 


<t>2 = + 2 ^3 “ ~ e(y-l) (t>2^/(3A0 2 ) 


( 6 . 8 ) 



and 
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• • 

- ~ < l , 2^^3 “ eW 2 

♦n_1 ■ -St+N-1 + 2C *N-2 " *N-l )(1 ' 6< ’ , " 1) *N-l )/(34e2) 

- 8£C *N-a * *N-2 K *N-1 - *n- 2 )/C9462) - eW N-l * C6 ' 9) 


For the phenomenological model,. the burning rate 
function becomes 

= nw(<j> i+1 - <|> ) 2 /C4 a 6 2 ) 2 < i < N-l . (6.10) 

Combining (6.8) and (6.10) one obtains 

v*2 = 4nw(<f>3 - <j>2 ) 2 /(9A0 2 ) 

W N-1 = 4n H.^N-l " jsj_2 ) 2 /(9A 02 ) (6.11) 


respectively . 

For the vaporization model, the burning rate 
function becomes 


w ± = w((l -e4>^/2)(l + (^^ + 2 

- <j>. . ) 2 /( 2 A 6 u t ) 2 ) lu - 1 )/e 2 2 < i < N-l (6.12) 

X” X i_j 


and, along the boundary, yields 
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w 2 = w((l - Jge* 2 ).(l + (2(^ - <J> 2 )/C3Aeu L )) 2 ) 3 ' 11 - l)/e 2 

w N _ 1 =w( d-^e4> N _ 1 ) Cl+C2((f. N _ 1 -({> N _ 2 )/(3Aeu L ) ) 2 ^-D/e 2 , (6.13) 

accordingly . 

Using either the phenomonological or the vaporization- 
limited combustion model, .(6. 5) and (6.9) can be solved by 
the Runge-Kutta method to obtain <|>^, i = 2,3>..«, N-l. 

Then the modal amplitudes can be obtained by using the 
Fourier -cosine transform 

IT 

f.(t) =2 / <j)COsi0de/Tr. (6.14) 

1 o 

These integrals must be computed numerically since <p is 
known only at the grid points. 

Some typical results for standing waves are presented 
in Table 3* They are intended to illustrate the 
accuracy of the two-term Galerkin solution and to compare 
the results associated with the phenomenological and 
vaporization-limited combustion laws. The column labeled 
PG indicates the results with a two-term Galerkin solutidn 
using the phenomenological model, the column labeled PF 
contains the results of a finite-difference solution of 
these equations, and the column labeled VF presents results 
of a finite difference solution using the vaporization- 


limited model 
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Tab le 3 


Stability Limits For Phenomenological and 
Vaporization-Limited Combustion Model 



PG 

PF 

VF 

e 

n 

n 

U L 

.i 

46 

28 

1.3 

.05 

91 

56 

2.1 

.01 

452 

286 

5.2 


It can be seen that the order of magnitude of the 
interaction index n required for instability for both the 
finite difference and Galerkin methods are roughly the same 
but the stability boundaries predicted by the Galerkin method 
are approximately twice as high as those predicted by the 
Finite-difference method. It is to be expected that the 
Galerkin method will lead to higher stability limits than 
the use of an exact solution procedure. This can be explained 
as follows. The instability mechanism is basically a feed- 
back process. Consider the form of (6.3) associated with 
standing waves in an annulus. This is 

M • • • 

f l + - f l + f l + 2e ( f i f 2 + f 2 f l^ + 2£ 51 nf i f 2 = 0 

f 2 + wf 2 + 4f 2 - efjfj - kevmf* = 0. (6.14) 
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For clarity the gas-dynamic nonlinearities will be neglected 
to give 


?! + wf 1 + f 1 + 2 ewnf 1 f 2 = 0 
f 2 + wf 2 + 4f 2 - ^ewnf* = 0. (6.15) 

A one-term Galerkin solution would lead to the 
equation 

f ^ wf + f ^ = .0 . 

This would indicate unconditional stability. Instability 
can arise only when the presence of the last term in (6.15b) 
causes f2 to grow and this then causes f^ to grow because 
of the last term in (6.15a). The growth of f2 also produces 
the growth of higher modes which in turn causes additional 
growth of f i . The contributions of these higher modes are 
neglected in a two-term Galerkin analysis and thus the energy 
input due to unsteady burning is underestimated. For this 
reason the Galerkin analysis will overestimate the stability 
of the system. This overestimation will decrease as the 
number of terms retained is increased. An example of this 
can be seen by comparing the one- and two-term analyses. 

A one-term solution predicts that the system is always 
stable. A two-term solution predicts the correct qualita- 
tive behavior but overestimates the system's stability by 
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roughly a factor of two. It is to be expected that further 
increases in accuracy can be obtained by keeping more terms 
in the assumed solution but that this will not seriously 
affect the quatitative prediction. 

In Chapter 5, it was found that equating dw/du 2 the 
phenomenological and vaporization-limited combustion laws 
at u 2 = 0 lead to the formula 

n = l/(4e 2 u£). (6.16) 

Assuming that the actual data can be fitted to a formula of 
the form 


n = C/(4 e 2 u 2 ) (6.17) 

the data give in Table 3 provide an opportunity to estimate 
C. The results are shown below. 


e 

.1 

.05 

.01 

c 

1.89 

2.47 

3.09 


It appears that a value of C = 2.5 would give accept- 
able accuracy over this range of e. 

For a given value of u^ the n computed from ( 6 . 16 ) 
will overestimate the burning rate. Thus it might be 
thought that C should be less than unity. The stability 
criterion was, however, based on the magnitude of the modal 
amplitudes. Inspection of Tables 4 and 5 shows that the 
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vaporization-limited combustion model involves the higher 
modes to a greater extent than does the phenomenological 
combustion model. Thus, a given value of Ul can be associated 
with lower individual modal amplitudes in the former case 
than in the latter. These two effects interact and calcula- 
tion shows that C is greater than unity in this range of e. 
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TABLE 4 

MODAL AMPLITUDE FOR PHENOHENLOGICAL MODEL 
IN ANNULAR COMBUSTOR 


• m ce • m 

T 

RODE 1 

MODE 2 

MODE 3 

MODE 4 MODE 5 

MODE 6 

MOOE 7 

MODE 3 

O.GO 

l.OGO 

0.000 

-0.000 

0.000 -0.000 

0.000 

-0. 000 

0. 00 0 

1.25 

0.3 37 

0. 029 

0. 002 

0.000 

0.000 

-0.000 

-0.000 

-0.000 

2.50 

-0.687 

-0.018 

-0.000 

0.000 

0.000 

-0-000 

0.000 

0.000 

3.75 

-0.700 

0.064 

-0.008 

0.001 -o.ooo 

0.000 

0. 000 

-0.000 

5. CO 

0. 2C9 

-0. 043 

0.009 

-0.001 -0.000 

0.000 

-0.000 

o.ooc 

6.25 

0.751 

0.040 

-0.006 

-0.003 -0.001 

0.000 

0-000 

0.000 

7.50 

0.237 

0.025 

0.000 

-0.002 -0.001 

-0.000 

-0.000 

-0.000 

8.75 

-0.535 

-0.060 

0.015 

0.004 -0.002 

-0.000 

0.000 

0.000 

10.00 

-0.523 

0.106 

-0-023 

0.009 -0.003 

0.001 

-0.000 

0.000 

11.25 

0.187 

-0.088 

0.024 

-0.000 -0.003 

0.002 

-0.000 

-0.000 

12.50 

0. 585 

0.062 

-0.009 

-0.011 -0.005 

-0.001 

0.001 

0.001 

13.75 

0.153 

0.011 

-0.010 

-0.010 -0.006 

-0.00 3 

-3.002 

-0.001 

15. CO 

-0.440 

-0.065 

0.032 

0.007 -0.007 

-0.001 

0.002 

0.000 

1 6.25 

-0. 390 

0 .112 

-0.043 

0.020 -0.010 

0.005 

-0. 002 

0-001 

1 7.50 

0.183 

-0.105 

0.030 

0.005 -0.010 

0.004 

0.001 

-0.002 

1 8.75 

0.463 

0.078 

-0.004 

-0.017 -0.013 

-0.004 

0.001 

0.002 

20. CO 

0.085 

-0.009 

-0.023 

-0.020 -0.014 

-0.008 

-0.004 

-0.002 

21.25 

-0. 374 

-0.050 

0-045 

0.004 -0.015 

0.001 

0.006 

-0.000 

22.50 

-0.287 

0.101 

-0.048 

0.028 -0.013 

0.011 

-0.007 

0.00 3 

23.75 

0.184 

-0.107 

0.025 

0.015 -0-018 

0.005 

0.004 

-0.005 

2 5. CO 

0. 371 

0.089 

0.009 

-0.016 -0.019 

-0.010 

0.001 

0. 00 5 

26.25 

0.030 

-0. 031 

-0.0 37 

-0. 030 -0.020 

-0.012 

-0.006 

-0. 002 

27.50 

-0. 327 

-0.027 

0.052 

-0.005 -0.021 

0.005 

0.010 

-0.002 

0.006 

28.75 

-0.206 

0.081 

-0-044 

0.030 -0.023 

0.017 

-o. on 

3 0. CO 

0. 1 91 

-0.102 

0.013 

0.028 -0.023 

0.003 

0.008 

-0.006 

31.25 

0. 3C2 

0.098 

0.026 

-0.009 -0.023 

-0.017 

-0.003 

0. 004 

32.50 
3 3.75 

--V.SH 

-0.054 

-0.051 

-0.038 -0.024 

-0.012 

-0.004 

0.001 

0. GOO 

0.056 

-0.019 -0.024 

0.012 

0.010 

-0.005 

35. CO 

-0.141 

0.053 

-0.034 

0.027 -0.024 

0.020 

-0.013 

0.008 

36. 25 

0.207 

-0.095 

-0.005 

0.040 -0.023 

-0.003 

0.012 

-0.005 

37.50 

0. 25 2 

' 0.108 

0.046 

0.004 -0.023 

-0.022 

-0.008 

0. 001 

38.75 

-0.067 

-0.079 

-0.064 

-0.042 -0.022 

-0.007 

0.002 

0.005 

40.00 

-0.284 

0.031 

0.056 

-0.036 -0.020 

0.019 

0.005 

-0.008 

41.25 

-C. 085 

0.033 

-0.019 

0.018 -0.019 

0.016 

-0. 010 

0-007 

42.50 

0.236 

-0.087 

-0.029 

0.052 -0.017 

-0.013 

0.013 

-0.001 

43.75 

0.218 

0.120 

0.070 

0.023 -0.014 

-0.022 

-0.014 

-0.006 

45. CO 

-0.122 

-0.109 

-0.076 

-0.039 -0.011 

0.006 

0-012 

0.010 

46.25 

47.50 

-0. 292 
-0.029 

0.070 

-0.001 

0.049 

0.005 

-0.055 -0.007 
-0.001 -0.003 

8:SS1 

-0.007 
0. 000 


48.75 

0. 290 

-0.073 

-0.063 

0.058 

0.003 

-0.026 

0.007 

0.009 

5 0. CO 

0.192 

0.135 

0.098 

0.053 

0.011 

-0 .008 

-0.014 

-0. 01 3 

51.25 

-0.207 

-0. 152 

-0.081 

-0. 020 

0.016 

0.026 

0.019 

0.006 

52.50 

-0.328 

0. 134 

0.022 

-0.073 

0.026 

0.015 

-0.022 

0.006 

5 3.75 

0.055 

-0.060 

0.055 

-0.044 

0.034 

-0.028 

0. 022 

-0.015 

55. CO 

0.4 04 

-0.033 

-0.119 

0.037 

0.046 

-0.022 

-0.020 

0.01 2 

56.25 

0.151 

0. 140 

0.125 

0.094 

0.061 

0.038 

0.019 

0.008 

57.50 

-0.403 

-0.219 

-0.045 

0.054 

0.065 

0.025 

-0.013 

-0. 022 

58.75 

-0.412 

0. 276 

-0.092 

-0.036 

0.071 

-0.051 

0.C13 

0.012 

60.00 

0. 324 

-0.246 

0.183 

-0.117 

0.060 

-0.017 

-0.011 

0.022 
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TA8LE 5 

MODAL AMPLITUDE FOR VAPORIZATION MODEL 
IN ANNULAR COMBUSTOR 


T 

MODE 1 MODE 2 

MODE 3 

MODE 4 

MODE 5 

MODE 6 

MODE 7 

MODE 8 

0.00 

1.000 

0.000 

-0.000 

0.000 

-0-000 

0.000 

-0. 000 

0. 00 0 

1.Z5 

0.384 -0.193 

0. 0 39 

-0.007 

0.002 

-0.00 1 

0.000 

-0. 000 

2.50 

-0.455 

0.095 

-0.033 

0.009 

-0.002 

0.000 

0.000 

-o.ooo 

3.75 

-0.514 -0.111 

0.042 

0.018 

-0.007 

-0.003 

0. 002 

0. 000 

5. CO 

0. 060 -0.035 

-0.036 

-0 .029 

-0.020 

-0.012 

-0.007 

-0.003 

6.25 

0.443 

0.106 

-0.022 

-0.041 

-0.021 

-0.002 

0.006 

0.006 

7.50 

0.293 -0.225 

0.077 

0.020 

-0.043 

0.019 

0.007 

-0.013 

8.75 

-0. 324 

0.160 

-0.104 

0.071 

-0.049 

0 .034 

-0.024 

0.017 

10. CO 

-0.717 -0.029 

0.192 

-0.037 

-0.072 

0.042 

0.020 

-0.032 

1 1.25 

-0.224 -0.455 

-0.330 

-0.172 

-0.040 

0.041 

0.065 

0.047 

1 2.50 

1.281 

0. 558 

0.300 

0.131 

0.040 

-0.003 

-0.021 

-0. 030 

13.75 

2.176 -0.505 

-0.474 

0. 087 

0.172 

0.024 

-0.078 

-0.053 

15.00 

0.264 -0.967 

0.644 

-0. 327 

0.068 

0.076 

-0.105 

0.057 

16.25 

-2.9C7 

0.857 

-0.311 

0.100 

-0.053 

0.005 

0.019 

-0.031 

1 7.50 

-3. 416 -1.125 

0.382 

0.037 

-0.193 

0.052 

0.072 

-0.074 

18.75 

0.448 -0.794 

-0.551 

-0.324 

-0.057 

0.085 

0.110 

0. 056 

20.00 

3.619 

0.546 

0. 200 

0.077 

0.001 

- 0.024 

-0.039 

-0.044 

21.25 

2.747 -1.259 

-0.224 

0.106 

0.166 

-0-024 

-0.087 

-0.025 

22.50 

-1.442 -0.328 

0.340 

-0.271 

0. 145 

-0.042 

-0. 032 

0.063 

2 3.75 

- 3. 229 

0.136 

-0.113 

-0.068 

0.068 

- 0.080 

0.058 

-0.038 

25.00 

-1.329 -1.066 

0.005 

0.181 

-0.028 

-0.10 2 

0.012 

0. 059 

26.25 

1.794 

0. 024 

-0. 009 

-0. 08 5 

*0 • 086 

-0.083 

-0.067 

-0.049 

27.50 

2.388 -0.063 

-0.093 

-0.167 

-0.089 

-0.030 

0.022 
0. 069 

0.041 

28.75 

0.463 -0.881 

0-235 

0.090 

-0.114 

-0.006 

-0.027 

3 0. CO 

-1. 8C7 

0.272 

-0.195 

0.095 

-0.067 

0 .043 

-0.031 

0.021 

31.25 

-2.088 -0.294 

0.261 

-0.143 

-0.024 

0.075 

-0.052 

-0. 00 3 

32.50 

-0.063 -0.796 

-0.404 

-0.100 

0.077 

0.093 

0.025 

-0.040 

33.75 

2.115 

0.417 

0.250 

0.120 

0.079 

0.049 

0.033 

0.021 

35.00 

2.172 -0.647 

-0.302 

-0.018 

0.128 

0.062 

-0-032 

-0.059 

36. 25 

-0.447 -0.603 

0.422 

-0.241 

0.058 

0.049 

-0.079 

0.050 

37.50 

-2.586 

0. 4 34 

-0.193 

0.046 

-0.004 

-0.025 

0.035 

-0. 039 

38.75 

-2.044 -1.000 

0.179 

0. 126 

-0.127 

-0.044 

0.074 

O.OOC 

40.00 

1.142 -0.253 

-0.260 

-0.241 

-0.154 

-0.072 

-0.004 

0.038 

41.25 

2.348 

0.253 

0.055 

-0.083 

-0.089 

-0.061 

-0. 054 

-0.027 

4 2.50 

1.497 -1.140 

0.032 

0.196 

0.007 

-0.109 

0.006 

0.063 

43.75 

-1.721 

0.082 

0.011 

-0.091 

0.095 

-0.093 

0.077 

-0. 057 

45. CO 

-2.759 -0.014 

0.116 

-0.188 

0.105 

-0-038 

-0.020 

0.044 

4 65 25 

-0.807 -1.084 

-0.251 

0.126 

0.118 

-0.0 3 3 

-0.081 

-0.012 

47.50 

2.050 

0. 301 

0.193 

0.083 

0.050 

0.024 

0.010 

-0.001 

48.75 

2.514 -0.289 

-0.260 

-0.186 

-0.011 

0.062 

0.068 

0.023 

50.00 

0.206 -0.932 

0.405 

-0.043 

-0-117 

0.082 

0.014 

-0. 063 

5 1.25 

-2.255 

0. 397 

-0.263 

0.134 

-0.096 

0.066 

-0.051 

0.040 

52.50 

-2.287 -0.575 

0.312 

-0.076 

-0.106 

0.084 

-0.002 

-0.053 

53.75 

0.331 -0.709 

-0.439 

-0.200 

-0.000 

0.08 3 

0. 074 

0.017 

55.00 

2.477 

0.407 

0.229 

0.082 

0.038 

0.004 

-0.011 

-0.021 

56.25 

2.025 -0.858 

-0.236 

0.073 

0.141 

-0.002 

-0.072 

-0.031 

57.50 

-0.915 -0.377 

0. 322 

-0.252 

0.129 

-0.030 

-0.036 

0.060 

58.75 

-2.669 

0. 303 

-0.109 

-0.033 

0.057 

-0.068 

0.057 

-0.042 

60. CO 

-1.583 -1.060 

0.038 

0.179 

-0.051 

-0.099 

0.028 

0.055 
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CONCLUSIONS 

The primary objective of this study was the 
development of a new analytical technique to be used in the 
solution of nonlinear velocity-sensitive combustion 
instability problems. Such a method should be relatively 
easy to apply and should require relatively little 
computation time. 

In an attempt to achieve this aim, the orthogonal 
collocation method was investigated first. However, it was 
found that the results were heavily dependent on the location 
of the collocation points and characteristics of the 
equations. Therefore, the method was rejected as unreliable. 

Next, the Galerkin method, which has proved to be 
very successful in analysis of the pressure sensitive 
combustion instability, was considered. This method proved 
to work very well. It was found that the pressure waveforms 
exhibit a strong second harmonic distortion and a variety 
of behaviors are possible depending on the nature of the 
combustion process and the parametric values involved. 

Finally, a one-dimensional model provided further 
insight into the problem by allowing a comparison of 
Galerkin solutions with more exact finite-difference 
computations . 
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Some major conclusions of this research are as follows. 
(1) The form of (2.33a) is somewhat Insensitive to the 
specific assumptions used to derive it. This is demonstrated 
by the fact that Powell (1) developed an equation of a very 
similar form using a different set of assumptions. (2) The 
orthogonal collocation method is unsuited to solution of 
problems of the type under discussion here. (3) The Galerkin 
method is well suited to the solution of such problems. (4) 
Stability boundaries and pressure wave forms appear to be 
highly dependent on the parameters of the problem (interac- 
tion index, time delay, steady-state burning rate, etc.). 
Furthermore , there appears to be no clear pattern to the 
computed results. (5) The phenomenological burning-rate law 
employed in the majority of the work discussed predicts 
results which are qualitatively similar to those associated 
with the vaporization-limited burning-rate law used by 
previous investigators. The phenomenological law, further- 
more, can be used in. conjunction with the Galerkin method 
while the vaporization-limited law cannot. (6) The computer 
programs developed in the course of this work can be 
employed to determine stability boundaries and pressure 
waveforms without the expenditure of excessive computer 
time. This is important because the lack of clear trends 
discussed above make an individual analysis of each 
situation desirable. 



REFERENCES 


1. Powell, E. A., "Nonlinear Combustion Instability in 

Liquid Propellant Rocket Engines," Ph. D. Dissertation, 
Georgia Institute of Technology, 1970* 

2. Maslen, S. H. , and Moore, F. K. , "On strong Transverse 

Waves without Shocks in a Circular Cylinder," 

Journal of Aeronautical Sciences, Vol. 23 (1956), 

PP. 533-593. 

3* Priem, R. J., and Guentert, D. C., "Combustion Instabili- 
ty Limits Determined by a Nonlinear Theory and a 
One-Dimensional Model," TN-l409> NASA Lewis Research 
Center (1962). 

4. Powell, E. A., and Zinn, B. T. , "A Single Mode Approx- 

imation in The Solution of Nonlinear Combustion 
Instability Problems," Combustion Science and 
Technology, Vol. 3 (1971), pp. 121-132. 

5. Lores, M. E., and Zinn, B. T. , "Nonlinear Longitudinal 

Combustion Instability in Rocket Motors," Combustion 
Science and Technology j Vol. 7 (1973), pp. 245-256. 

6. Peddieson, J., Jr., Ventrice, M. B. , and Purdy, K. R., 

"Effects of Nonlinear Combustion on the Stability of 
a Liquid Propellant System," Proceedings of 14th 
JANNAF Combustion Meeting (CPIA Publication 292) , 
pp. 525-53# (1977). 

7. McDonald, G. H. , "Stability Analysis of A Liquid Fuel 

Annular Combustion Chamber," M. S. Thesis, TENN. 

TECH. UNIV., August 1979- 


56 



57 



TRAVELING WAVE 



Figure 1. Pressure functions for standing and traveling 
waves 



























































































APPENDIX A 


DERIVATION OF THE COEFFICIENTS 

In this section, the coefficients b . . . C , and d. . 

ij ij ij 

appearing in chapter 3 and 5 are derived and their numerical 
values are presented in Table 6, 7>and 8 respectively. 

In chapter 3 , for initial condition (3.7 ) , one assumes 

oo 

9 tt ^2 “ v 2 ^ 2 = m ^ 0 K m ^ r ) cos ) sin C2S 11 t) • 

By comparing coefficient yields 

H q (r ) = ^S 11 C2S| 1 J2 - Cy-DSJjJJ - * Sl iWr + ^g/r 2 ) 
H x Cr) = 0 

H 2 (r) = ^S 11 (2S2 1 J§ - (y-DS^Jf - 4S 11 J Q J 1 /r) 

H (r) = 0 for m = 3, 4. ...,<*». 

m 

Then expanding <J > 2 in a Fourier-Bessel series and expanding 

H (r) in a Bessel series leads to 
m 

b ij . = 2S2 j /^H 1 (r)J i (S iJ r)rdr/{(S 2 j - i 2 ) J^CS^ ) > . 
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This integral must be computed numerically. For the sake 
of brevity only the first five coefficients in each series 
are calculated and presented in Table 6. 

Table 6 

BESSEL SERIES COEFFICIENTS FOR ANALYICAL SOLUTION 


n 

b 0n 

b 2n 

0 

1. 11738-0. 37246y 

• 

1 

0.50462+0. 37912y 

0. 27181-1. 09482y 

2 

-0. 09246-0. 00802y 

-0. 11390-0. OIHOy 

3 

0. 05019+0. 00184 y 

0. 05479+0. 00208y 

4 

-0. 03290-0. 00067y 

-0. 03451-0. 00071Y 

5 

0. 02375+0. 00030 y 

0. 02451+0. 00033y 


In Chapter 5, equation (5.8) can be written as 


$ ^*0^0 ^*1^1 ^ 2^2 ®1^*3 ^ 4 (A-l) 

where <J>q = Jq, <{>i = Jocose, 4*2 = JgCOsSe 

<J>2 = J 1 sin0, ^ = J 2 sin20. 

By substituting (A-l) into the governing equation (5*6), 
yields 


R = D(4) 


CA-2) 
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where D is the nonlinear differential operator of equation 
(5.6). If CA-1) represented the exact solution to (5.6), r 
would vanish. Since this is not the case R has a finite 
value. The Galerkin procedure consists of making R orthog- 
onal to each of the <j>^'s. This leads to the equations 

f / 21T R<(. rdrd = 0 • i = 0 , 1 ,..., 4 CA-3) 

o o i 

and a set of five second order differential equations are 
generated. The angular part of integration can be performed 
in closed form while the radial part of integration must be 
computed numerically and leads to the integral 



°i = Voh Cr5J k Cs ki rirdr 

CA-4) 

where 

6 k - 2s kl /( ^- S kl " k2 ' J k^ S kl^ , ‘ 

CA-5) 


The numerical values of are listed in Table 7 for 
Y = 1.2. Also listed in this table are the functional forms 
of F i , s. In writing these the notation 

R i = + s n J i /r " l2j ±/ r2 > 1 = 0 > 1,2 (A-6) 

is used. Similarly the numerical values of 

■ e i /J 0 G ij J i Cs ii r)rdr 


(A-7) 



and 


fun =tionai tom 

orms or a 

iJ S llst 

° ln Table 8, 
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Table 7 


COEFFICIENTS APPEARING IN EQUATIONS (5-9) Y 


4.137 

1.042 


- 0.208 


R 0 J o + 2S§ It Jj2 

RjJj/2 + Sf 3 J] 2 + J 2 /r 2 

R 2 J 2 /2 + S§ 1( J£ 2 + 4j|/r 2 


-1.939 R 0 J 3 + 2S 0J S 1I J«Jj 

-2.312 RiJ o + 2s oi s n J o J i 

1.719 R 2 J 2 /2 + S 1]L S 21 J;J» + PJ^/r 2 

1.483 R 2 J 1 /2 + S 11 S 21 J]J» + 2J 3 J 2 /r2 

-2.785 R 0 J 2 + 2S 0 jS 23 J g J 2 

-3.039 R 2 J 0 + 2S 01 S 21 J^ 

1.132 R jJj/ 2 + S 2 lt Jj 2 - J 2 /r 2 

2.586 . S§ 2 JJ 2 

0.480 ^(Sh-Ji 2 + Jf/r 2 ) 

-0.196 S 21 J 2 2 /2 + 2J |/r 2 

-4.849 2S 01 S 1:i J*j; 

1.503 + 21^/r 2 

-0.576 ^(S 2 J Jj 2 - J 2 /r 2 ) 


-3-481 


2S nl S J«J: 
01 21 0 2 






BIS] 


Table 8 


COEFFICIENTS APPEARING IN EQUATION (5-12) 


wSSM 


1.000 

J 0 CSoir) 

1.293 


0.240 


-0.098 

JsCsfij * 2 / 2 

-0.176 

OslO 

ho 

l 

0.061 

-Vf 

0.04g 


1.000 

J 1 CSur ) 

-1.212 

s oi^u J o J i 

0.927 

^5(8^2^21'^ i' 

0.165 

“ J 0 J 1 

-0.199 

“^2 J 2 J 2 

1.000 

J 2 (S 2 ir ) 

-1.741 

s oi s 2 i j ; j ^ 

-0.223 

^ s h j i 2 - 

0.237 

-J 0 J 2 

-0.175 









APPENDIX B 

PROGRAM FOR STABILITY BOUNDARY CURVES 
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SSET LIST FREE 
C 

C KIN-WING WONG JUN 30 1978 

C 

C CCMBUSTION INST ABILITY 

C 

C. METHOD GALERKIN METHOO 

C« ThIS PROGRAM GENERATE THE STABILITY CURVE. 

C 

C INPUT INFORMATION 

C 

C THE PROGRAM USED THE FORMAT FREE INPUT OUTPUT. A 


c 

COMMA SERVES 

AS A DELIMITER ANO THE INPUT VARIABLES 

c 

ARE ENTER IN 

THE FOLLOWING ORDER 

c 

T I 

INITIAL TIME 

c 

TF 

FINAL TIME 

c 

DT 

STEP SIZE FOR TIME 

c 

MB 

BURNING RATE 

c 

r 

tdelay 

DELAY TIME 

U 

c 

El 

INITIAL EPS 

c 

EF 

FINAL EPS 

c 

OE 

STEP SIZE FOR EPS 

c 

OA 

ASSUME POSITION FOR N 

c 

F 0»FD 

ARE THE INITIAL CONDITION MODIFIER ANO 

c 

c 

c 

ANOTHER SETS 

HAVING THE VALUE EITHER 1 OR 0 
OF DATA CAN 8E ENTER BY ENTERING 

c 

r 

E I # EF»DE»I 

DA>FO AND FD. 

L 

Cl ME NS I ON DY<1O)#FO(5)#0F<5)#GI1O#5OO)*Y< 10) 
COMMON EP# RK»mB*EPANV 


COMMON /SLIST/S0»S1*S2 

COMMON /BLK1/C1»C2»C3»C4»C5»C6»C7#C8#C9»C10»C11»C12 
1* Cl 3* C14. C 15#C16»C 17 

LOGICAL PA SS»GRO h 
DATA MAXITR#EPS/20*.l/ 

DATA PI/3. 1*159/ 

READC5*/) TI#TF#OT»WB*TDELAY 
PFINT*//»TI*TF»DT*W8*TDELAV 
UELAY=PI*TOELAY 
K=TDELAY/DT 
K 1= K •* 1 

IF CK1.GT.5C0) GO TO 105 
If (K.LT.l.AND.NCF.Efi.l ) R=1 
NSTEP = (TF“ TI )/D T 
N 1 1 0 

HhALF=DT*.5 
N FRO 8 = 0 

REA D{ 5#/»£ND=9S9 ) EI»EF*DE#DA»FO*CF 
NFROB=NPRO04l 


9 



O O O 
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F FI NT *// .NPFQB 
PFINT *//#F0*DF 
HSTEP*(EF-EI)/0£ 

Ef =E I 

IF (EP.LT.O) EP«.Ol 
AN¥=DA 

FASS=. FALSE. 

CC IOC ISTEP =0#MSTEP ' 

CC 80 ITR*l*MAXITR 
EFAN V=EP*AN V 

IF (EPANV-EC.O) EPANV=ANV 
CC 10 I = 1»K 
CC 10 J-l.N 
1C G(J*I)=0. 

CC 20 I*WN 

2 C Y(I)=G. 

INITIAL CONDITION 

YC2>=F0tl>m* > = F0C2>;Y( 6>=F0(3)M(8)-F0< 4 >;y< 10J=FO 
1(5) 

Y (1 )*CFC1 )*SOA Y( 3)=DF<2)*SHY<5>=0F(3)*$2 
Y(7) a 0F(4)*SliYCS) = DF(5)*S2 
CALL DLYFUN <N*Y*6»K*EPS#Kl) 

T = T I 
C 

CC 40 I *1 * NSTEP 

CALL RKDELY(N*T #Y»DY*DT#G*K#EP»411Q»HHALF*K1 ) 

CC 30 J=2.N*2 

IFC ABSCYCJ )).GT.2) GO TO 60 

3 C CCNTINLE 

4 C CCNTINUE 

GFOW=. FALSE. 

IF (PASS) GO TO 50 
AN VD MN = AN V 
ANV= ANY+OA 
GC TO 80 
50 CCNTINUE 

A F VD WN = ANY 
GC TO 70 
60 CCNTINUE 

P ASS=. TRUE. 

GFON-.TRUE. 

AFVUP-ANV 

7 C CCNTINLE 

IF ( A8S(ANYUP-ANYDNN).LT.EPS) GO TO 90 
ANV= ( ANVUP + ANVD nN)*«5 

8 C CCNTINUE 
90 CCNTINUE 

NNU = AN V 

PFINT *//»EP*NNU#ITR 



ooo O OOOOOOOO 
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IF ( ITR.EQ.MAXITR ) PRINT ✓ /#•••• WARNING *** THE 
1 ANSWERMAY NOT CONVERGE* 

EF=EP«0E 

ANV=ANVUP/2#ANVDWN=0 
ICC CCNTINUE 
GC TO 9 

1C5 PRINT //* * DELAY TIME TOO LARGE • #TDEL AY 
GC TO 999 

11C FFINT *//* * ILLEGAL DELAY FUNCTION* 

999 S10P 
END 

SUBROUTINE DIFFUN CT* Y*D Y ) 

IHS SECTION PROVIDE A SET OF FIVE SECOND ORDER 
EOLATIONS 
WHERE 

F0= YC2 ) G1 = YC8) 

F 1 = YC 4 ) G2=VC1C> 

F 2 3 Y( 6) 

DIMENSION YC10)#CYC10> 

C C MM ON EP* FIC*W8#EPANV 
CGMMON /SLIST/S0#S1*S2 

COMMON /BLKl/ C 1#C2#C 3*C4»C5#C6#C7 #C8#C9# Cl 0*C 1 1*C 12 
1* Cl 3* C14# C15»C16#C1? 

CYC2)=YC1) 

DYC 1) = -S0*SC*YC2)-EP*CC1*YC2) *YO)*C2*CYC4)*YC3 )*YC8) 
1* YC 7 ) ) ♦C3*CYC10>*YC9)*YC6>*YC5)>> 

CYC4 )= YC 3) 

DY( 3) 3 ”S1*S1*YCA)“EP*<C4*Y<2)*Y(3)*C5*YCA)*Y<1) ♦ 

1 C 6 * C Y C 8 )* YC9J* YC4 )*YC5 ) >*C7*CYC6 )* Y( 3 )♦ Y < 1 0 >* YC 7 ) ) 1 
CYC6)=YC5) 

CYC5 )a-S2*S2*YC6>-£P*CC8*YC2>*YC5)*C9*Y<6 )*YCl >♦ CIO* 
1CYC8)*YC 7) * YC 4 ) * YC 3)) ) 

C Y C8 }- Y C7 ) 

C Y (7 >»-Sl*Sl*YC8 >-EP*CC4*YC2>*YC7)*C5*YC8 J*YC1 >♦ C6* 

1C YC4)*YC9)-YC8)*YC5))*C7*C YCIO )*Y( 3 >-VC6 > *YC 7 )) ) 

DYC 10)=YC9) 

CYC9 )=-S2*S2*YC10)-EP*CC8*YC2)*YC9>*C9*Y< 10>*YC1) CIO 
1*CYC8)*YC3)*YC4>*YC7>>> 

IF (RK.EQ. 2) RETURN 

THIS PROVIDE THE COMBUSTION TERMS 

CYC l)=DYC 1 )-WB*< YC1 )*EP ANV* CC 1 1 *Y C 2 )* YC 2 >*C 12*C YC 4 ) *Y 
ICO* YC8 ) *Y C8 ) )* C13*CYC10 )* Y Cl G > ♦ Y C 6 >* Y C 6 > > > ) 

CYC 3 ) = DYC3)-WB*CY( 3 >*EPA NV *CC 14 *YC 2 >* YC 4 )*C1 5* C YC4 )*Y 
ICE)* YC 8 ) * YC 1 0 ) ) ) ) 

CYC 5 ) = DYC5 )-W8*C YC5 )* EP AN V* (Cl 6* C Y C 4 )* YC 4 )- Y C 8 ) *Y C 8) ) 

1 *C17*YI2 )*YC6>> ) 



etc 7)=DY(7 )-M8*(Y(7)4EPANV*(C14*Y<2 )* Y<8 )*C15*( Y< 4 )*Y 
1(10) Yl8)*YC6) >)) 

0Y(9)=DYC9 )-WB*( Y< 9 )4 EP AN V *(C 1?*Y(2)*Y( 10)42*C16*Y(4) 
1*Y(8))) 

RETURN 

END 

SIBROLTINE RKDELY(N*T#Y*DY#H*G*K»EPS* *# HHALF *K 1 ) 

THIS SECTION PERFORM A FOURTH ORDER RUNGE'KUTTA METHOD 
KITH TIME DELAY FUNCTION 

DIMENSION YC 10 ) *0 V(10 )#Y2C10)»Y3( 10 )»G( 10*500) 

CCMMON EP# RK * MB *EP AN V 

COMMON /8LKI/ C 1»C2*C 3*C4*C5* C6*C7 *C8* C9* Cl 0*C 1 1 *C 12 
1 # Cl 3# C14* C15#C16*C17 

CALL OIFFUN(T*Y*OY ) 

CC 10 I =1 * N 

Y2 ( I )=Y (I ) ♦ HHALF*<DY(I) 4 G(I*1)> 

G(I*1)=(G< 1*1)46(1*2))/ 2. 

CONTINUE 
T = T 4HH ALF 

CALL C I FFUN ( T* Y 2 #D Y ) 

00 20 I = 1*N 

Y3CI)=Y(!) ♦ HH ALF * (D Y< I ) 4 G(I*l)) 

Y 2 ( I ) = Y2CI ) * 2»* Y 3< I ) 

CALL OIFFUNC T*Y3*DY) 

CC 40 1=1* N 

Y 3( I )= Y < I ) ♦ H* ( D Y ( I ) ♦ G(I*1)) 

IF(K.LT.l) GO TO 40 

DC 30 J=l* K 
6(1* J) = G(I # J4l) 

Y2( I ) = Y2( 1 ) 4 Y 3( I ) 

T = T 4 HH ALF 

CALL 0 1 FFliN (T *Y3*DY) 

CC 50 I = 1 * N 

YCI) = (Y2CI) - Y ( I ) 4 HHALF* (OY(I ) ♦ G(I*l)))/3. 
CONTINUE 

ENTRY DLYFU'N ( N * Y *G * K# EPS#K1 ) 

TUS GIVES THE FORM OF DELAY FUNCTION G 

IF (K.LT.l) RETURN 
IF (RIUEQ.2) RETURN 
G( 1 *K 1 )= MB* 

1(4)4 Y( 6 ) *Y(8 ) )4 
G ( 3 #K l ) = wfi* 

1(€)4 Y(8)*Y(10))) 

G(5*K1 ) = M8*EPANV*(C16*( Y( 4 ) *Y (4 )-t ( 8) * Y< 8 ) )4 C17 * Y( 2 )* Y 
1 ( 0 ) 

G (7 * K 1 )= MB* EPANV*(C14*Y(2)*Y(8)4C15*(Y(4)*Y 

1(10) Y( 8) *Y ( 6 ) ) ) 


EPANV*<C11*Y(2)*Y(2)4C12*( Y(4 )*Y 
Cl 3*( Y(1 C )* Y< 10)4Y(6 )*Y(6) )) 
EPANV*<C14*Y(2 )* Y(4 )4C15*( Y(4)*Y 



G (9 » K1 )=WB*EPANV*(C17*Y(2)«Y( 1Q)*2*C16*Y< 4>*YC€ ) ) 

RETURN 

END 

BLOCK OATA 

COMMON /BLK1/C l#C2*C3*C4*C5*C6*C7*C8*C9»C10*Cil#C12 
1»C13#C14, C15#C16*C17 

C EMM ON /SL IST/S0?S1*S2 
CAT A SO *SWS2/ 3. 83 171 >1.8*1 18 *3.0S424A 
CAT A Cl#C2*C3#C4/4.l373*l.0423#-.2084#-i. 9 394/ 

CATA C5»C6»C7*C8#C9 /-2. 3123# 1-7187 *1 .4828#-2-785# 
1*3.0368/ 

CATA C10»C11»C12*C13/1.1318*2.586#.48C*".196/ 

DATA C 14*C l 5*C16#C17/“2. 424 3* 1.853444 5* -.447 *-3.481/ 
END 
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PROGRAM FOR MODAL AMPLITUDE AND PRESSURE 
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£ SET 
S3 IND 
FILE 

C 

c 

c 

c 

c 

c 

r 

C 

C 

C 

C 

C 

c 

c 


LI STUTO0IND 

= FROM *F0STLI8/= FREE 

1(KInC=DISK*MAXR£CSIZE=15*8L0CKSI2E=420*AREAS=4 
1*AREASI2E = 1400*FILETYPE = 7 ) 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


KIN-WING WONG 


COMBUSTION INSTABILITY 


JUN 30 1978 


METHOD: THIS PROGRAM USED THE GALERKIN METHOO TO 

GENERATE FIVE SECOND ORDER DIFFERENTIAL EQUATIONS. 

A FOURTH ORDER RUNGE-KUTTA METHOD IS USED TO SOLVE 
THE RESULTING EQUATIONS. 

THE PROGRAM GIVES THE MODAL AMPLITUDE ANO THE PERTUBED 
PRESSURE IN GRAPHICAL FORM. 

THERE ARE TWO TYPES OF INPUT DATAS. THE FIRST SET OF 
INPUT DATA USED FORMAT FREE INPUT. A COMMA SERVES 
AS A DELIMITER AND THE DATA ENTER AS FOLLOw: 

. T I INITIAL TIME 
TF FINAL TIME 
DT STEP SIZE FOR TIME 
NP PRINT FREQUENCY 

NO F l FOP THE DELAY TIME APPROCHING 0 
0 OTHERWISE 

IGAS l FOR GAS-DYNAMIC TERMS 
0 OTHERWISE 

FO * INITIAL CONDITION MODIFIER 
DF El THER 1 OR 0 

R RADIUS OF THE CYLINDER 

THE SECOND SET OF DATA USED THE NAMELIST (LIST) 

T OPTION. THE DATA CAN 8E ENTERED IN 
HOWEVER THE FOLLOWING DATA MUST BE 


c 

INPUT OUT 

c 

ANY ORDER 

c 

PROVIDED 

c 

TDELAY 

r 

w* 

EP 

c 

wB 

c 

ANV 


DELAY TIME IN MULTIPLE OF PI 
ORDER PARMETtR f 
STEADY BURNING RATE 
INTERACTION PARMETER 
PLOT T FOR PRESSURE PLOT 
F OTHERWISE 

IF PLOT IS TRUE* ONE MORE DATA (NUMPT) IS NEEOED TO 
PROVIDE THE NUMBER OF POINT WITHIN 0 AND 2*PI. 

C1MENS10N XXC100)*YY(100)*DY(10)*P(5)*FO(5)*DF(5)*G 
1(10* 5C 0 ) * Y ( 10) 

CCMMON EP* IGAS* RK*W3*EPANV 
CCMMON /SL1ST/SC*S1*S2 

COMMON /3L K 1/ C l *C 2 *C 3 * C 4 *C 5* C6 *C7 * C8 *C9 * Cl 0* Cl 1* C 1 2 
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WC13*C14, C15>C16*C17 

COMMON /BLK2/C02»C03*C04»C05»C06*C07*C12P#C1 3P»C1 4p 
1 > Cl 5P » C22>C23*C24*C25»GAMA 

COMMON /BL K 3/ 3 JO# B J1 # B J2 # COS Q#C 0S2 Q* SINQ# SIN 20 
DATA PI/3.14159/ 

LOGICAL PLOT 

NAMELIST /LIST/ TDEL AY # E P# W8* AN V » PL OT 
READ (5./) TI#TF»DT#NP»NOF*IGAS*FO#QF#R 
PRINT *//#TI#TF*OT#NP*NDF 

IF CIGAS.EO.O) PRINT //**GAS CYNAMIC TERM IS OFF* 

PRINT *//#FO#DF#R 

EFS=l.£-6 

CALL BESJC SO*R# 0 # B J 0* E PS * IER ) 

CALL 3ESJ(S1*R*1#BJ1»EPS#IER) 

CALL BESJC S2*R#2»BJ2#EPS»IER) 

999 READ(5*LIST#END=99) 

REHINO 1 
TCLY=PI*TDELAY 
EFANY = EP*ANV 

IF (EPANV.EO.O) EPA NV = ANV 
WR I T E ( 6*LI ST ) 

N = 1 0 

riHALF = DT*« 5 
K=TDL Y/DT 

K 1 =K 4 1 

IF ( K 1 .GT. 500) GO TO 500 
IF (K.LT.l.AND.NDF.EO.l) K=l 
OC 30 I = 1* K 
CC 30 J=1 » N 
30 S(J#I>=0. 

DO 10 I = 1*N 
10 Y C I ) = 0 • 

C INITIAL CONDITION 

C 

T =T I 

Y(2)=F0C1);y( 4)=F0C2)AY(6)=F0C3)J Y(3)=F0T4)i Y(I0)=F0 
1(5) 

YC1 ) =DF(1) *SO;Y( 3)=DF< 2)*S1 J Y ( 5)=0F ( 3 )*S2 

Y(7)=DFC4) *Si; Y(9)=DF(5)*S2 

NSTEP=(TF-TI )/DT 

IF ( .NOT. PLOT) GO TO 40 

CALL PRESUE ( Y*P ) 

WRITE(l) T*P 
4 C CCNTINUE 

PRINT //^"FUNCTION" 

*RITE(6#100) T*( Y(I )»I=2#N#2) 

45 CONTINUE 

CALL OLYFUN ( N » Y *G #K » EPS » K1 ) 

DC 20 1 = 1 » NSTEP 

CALL RKDELY(N»T#Y»DY#DT#G#K#£P* &400#HHALF,K1 > 

DO 55 J=2#N#2 
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IF(ABS(Y(J)).LT.2> GO TO 55 

R H = 2 

55 CONTINUE 

IF (MOD(I>NP).NE.O) GO TO 20 
IF (.NOT. PLOT) GO TO 60 
CALL PRESUE( Y*P ) 

WRITEQ) T » P 

60 *RI TE(6#100) T# (YCJ)* J=2*N»2 ) 

65 CONTINUE 
2 C CONTINUE 

IF (RK.GT.O) PRINT //*• UNSTABLE* 

IF (.NOT. PLOT) GO TO 70 
ENDFILE 1 
REWIND I 

PRINT //,' PRESSURE « 

READ(5*/) NUMPT 
D6=PI *2./NUMPT 
NLM=NUMPT4 1 

temp = c;xx( n=o; 

CC 75 I = 2» NUM 
TEMP = TEMP4DQ 
75 XX(I)=TEMP 

98 R E A 0 ( 1 #ENO = 70 ) T,P 
PRINT //»"TIH£= " * T 
SUM Y = C • 

DC 80 1 = 1# NUM 
QC=XX( I ) 

Y Y< I ) = *GAMA*(P( 1)*P(2)*COS(QO)4P(3)*COS(2*QQ)4P(4) 

1 * SI N ( 3 Q )4 P(5)*SlN(2*flO) ) 

SUM Y = S UMY* Y Y( I ) 

60 CONTINUE 

IF (SUMY.NE.O.) CALL PL OT 20C X X # Y Y#NUM * 1 00 # 1 00) 

GO TO 98 
70 CONTINUE 
RMO 

GC TO 999 

AOO PRINT *//, 'ILLEGAL OELAY FUNCTION* 

GC TO 99 

5 OC PRINT *//» • DEL A Y TIME TOO LARGE* 

99 STOP 

100 FORMAT (X#8E12.5) 

11C F0RMATC13X.7E12.5) 

ENO 

SUBROUTINE DIFFUN (T#Y#DY) 

C1MENSI0N Y (10 ) # OY ( 10 ) 

COMMON EP# I GAS * RK#W3#£PANV 
COMMON /SL IST/S0#S1»S2 

CGMMON /BL K 1/ C 1 #C2# C 3# C 4# C 5# C 6# C7 # C8 #C 9# C 1 0 #C 1 1 # C 12 
1# Cl 3# C 14# C15#C16#C17 

C 

CY(2 )=Y(1 ) 
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CYC 1 )=-S0*S0*Y(2 )-EP*IGAS*(Cl *Y(2)*Y( 1MC 2*(Y(4)*Y( 3) 
l4Y(3)*YC7) > ♦ C3*CY(10)*Y(9)4Y(6)*Y(5>>> 

DY( 4) = Y( 3) 

DY(3 >=-Sl*Sl*Y<4 )-EP*IGAS*(C4*Y<2)*Y( 3 >4 c 5*Y( 4 > *Y< 1 ) ♦ 
1 C6*(Y(81*YC9)*YC4 >• Y ( 5 ) )♦ C7*C Y( 6) *Y( 3 )♦ Y C 10>*Y 

2(7))) 

0 Y ( 6 )= Y(5 ) 

DY(5 ) =-S2*S2*Y< 6 >-E P* I GA S * ( C8 * Y( 2 ) * Y( 5)4 C 9*Y ( 6 > *Y( 1 >4 
1 C10*(Y(8)*YC7)-Y(4>*Y(3>)> 

DY(8)-Y(7) 

DYC 7> = -Sl*Sl*Y<8)-£P*I GAS*(C4 *Y< 2)* Y(7 >4C5*Y(8 )*YQ >4 
1 C6*(Y(4)*Y(9)-Y(8)*Y(5> HC7* (Y<10)*Y(3>- VC6 >*Y(7 ))) 
OYC 10 ) = YC9 ) 

0YC9)=-S2*S2*Y( 10)-EP*IGAS*(C8*Y<2)*Y(9)4C9*Y(10>*Y<1) 
1 C10*(Y(8 )*Y( 3 )4Y(4)*Y(7 )) ) 

IF (RK.EQ. 2) RETURN 

DYC1 ) = DYC1 )-H0*(Y(1)*EPANW*(C 11* Y( 2 ) * Y( 2 ) 4C1 2 *( Y( 4 ) *Y 
1(A)4 Y( 8 ) *YC8) )♦ C13*( Y(10 )* Y(10> 4Y(6 >*Y(6 > )) ) 

DY( 3 ) = DY( 3)-W8*( YC 3)*EPANtf *(C 14*YC 2) *Y( 4 )*C1 5*{ Y( 4 )*Y 
1(6)4 Y( 8 )*Y( 1 0 )) ) ) 

CYC 5 > = DY(5 >-WB* ( Y< 5 )♦ EP AN V* (C 16*CY( 4) *Y( 4 > "Y ( 8 > *Y( 8) > 

1 4C17*Y(2)*Y(6)>) 

DY(7)=DY(7 >-wB*C YC7 )4£P A N V* <C 1 4* Y( 2 ) * Y ( 8 ) *C 15 *( Y< 4 ) *Y 
1(10) Y(8)*Y(6)))> 

C Y(9 ) = 0Y(9 >-WB* ( YC9 )♦ EP AN V* (C 1 7* Y( 2 )* Y( l 0 )♦ 2*C 1 6*Y( 4) 
1*Y(8 ) ) ) 

RETURN 

END 

SUBROUTINE RKDELY<N*T»Y*DY*H»G*K*£PS* **HHALF*Kl) 
DIMENSION Y(10)#DYC10)*Y2(10) *Y3C10 )#G( 10*500) 

COMMON EP*IGAS> RK* N8»EPANV 

COMMON /8LK1/ C 1 ,C2 »C 3* C4,C 5* C6#C7 » C 8,C 9# Cl 0#C 1 1 » C 12 
1 » Cl 3 » C 1 4* C15*C16*C17 

CALL D I FFUN ( T# Y* DY ) 

DO 1 I = 1*N 

Y2(I) = Y(I) 4 HH ALF * (D Y( I ) ♦ G(I*1>> 

G(I*1)=(G( I»l)4G(I»2))/2 

1 CONTINUE 

T = HHHALF 

CALL DIFFUN(T*Y2*DY) 

DO 2 I = l#N 

Y3(I)=YCI) * HHALF*(DY(I) 4 G(I>1>) 

2 Y2 ( I )= Y2 ( I ) 4 2* Y 3< I ) 

CALL 0IFFUN( T * Y 3» DY > 

DO 3 I =1»N 

Y 3 ( I )=Y ( I ) 4 H*(DY(I) 4 6( I#1 ) ) 

IF(K.LT.l) GO TO 3 
DC 10 J = 1 * K 
10 G ( I » J ) = G( I » J *1 ) 

3 Y2(I )=Y 2( I ) 4 Y 3( I ) 

T=T4HHALF 



non 


in 


CALL DIFFUN CT ,Y3#0Y> 

DC 4 1 = 1, N 

Y ( I ) = (Y2(I) “ Y< I ) ♦ HHALF*(DY<I ) ♦ GCI#l>))/3. 

4 CONTINUE 

ENTRY DLYFUN ( N,Y*G,K,EPS,K1 ) 

IF (K.LT«1 ) RETURN 
IF (RK.E0.2) RETURN 

. G< 1 »K 1 )= W8* EPANV*(CU*Y(2)*Y(2)*C12*( Y(4)*Y 

1(A)* Y(8)*Y(8>>* C13*(Y(10)*Y(10)*Y(6)*Y(6))) 

G( 3»K 1)» MB* EPANV*(C14*Y<2 >*Y(4 )*C15*(Y(4)*Y 

1 C € > ♦ Y(8)*Y(10>>> 

G(5,K1)=W8*EPANV*CC16*(Y(4)*Y(4>-YC 8) *Y( 8 ) )♦ C17* Y ( 2 )* Y 
1 ( 6 )) 

G(7,K1)= MB* EPANV*(C14*Y(2)*YC8)*C15*< Y(4)*Y 

1(10) Y(8)*Y(6))) 

G(9,K1 )»WB*EPANV*(C1 7*Y(2)*Y( 10)*2*C16*Y< 4>*Y<8 >) 
RETURN 
END 

SUBROUTINE PRESUE(Y,P) 

THIS SECTION PERFORM THE PRESSURE CALCULATION 

DIMENSION Y ( 10 )> PCS ) 

COMMON EP,IGAS, RK,wB,EPANV 

COMMON /8LK2/C0 2,C03, C04,C05,C06,C07,C12P»C13P*C14P 
l»C15P, C22,C23,C24,C25»GAMA 

COMMON /BLK3/BJO,BJ1»BJ2#COSQ ,C0S2Q,SINQ,SIN2Q 
OFO = Y(1>;OF1 = YC 3>f DF2=Y(5)*DG1=Y(7);DG2=Y<9) 

F0=Y(2); F1=Y(4>* F2 = Y(6); G1 = Y(8)) G2=Y(10) 
P(1)=BJ0*( DF0*EP*CC02*F0*F0*CQ3*CF1*F1*G1*G1 >*C04* (F2 
1*F2*G2*G2) *C05*DF0*DF0*C06*CDF1*DF1*DG1*0G1 )*C07* 

2(CF2*0F2* OG 2 *DG2 ) ) ) 

P (2 > = BJ1*(DF1*EP*(C12P*F0*F1*C13P*CF1*F2*G1*G2 >*C14P 
1 *DFO *DF 1* C15P*(DF1*DF2*DG1*DG2>) ) 

P( 3)=8J2*(DF2+EP*(C22*F0*F2*C23*(F1*F1-GI*GI )*C24*0F0 
i*DF2* C25*(DF1*DF1-0G1*0G1 ))) 

P(4) = 8J1*(DG1*EP*(C12P*F0*G1*C13P*CF1*G2-GI*F2)*C14P 
1 *CFO *DG 1* C15P*(0F1*0G2“DF2*DG1 ) )) 

PC5)=6J2*( DG2*EP*CC22*F0*G2*2.*C23*F1*G1*C24*DF0*DG2* 

1 2.*C25*0F1*DG1 )) 

RETURN 
END 

BLOCK DATA 

COMMON /BL K 1 / C 1 ,C2,C 3, C 4# C 5, C 6, C7 , C8 ,C 9 , Cl 0 ,C 1 1 ,C 12 
1*C1 3*C 14, C 1 5 ,C 16,C 17 

COMMON /$LIST/S0,S1,S2 

COMMON /8LK2/C02»C03,C04,C05,C06,C07*C12P»C13P,C14P 
l # C 1 5P , C22,C23,C24,CZ5»GAMA 

DATA SO, SI ,32/3.83171,1.84113,3.05424/ 

DATA Cl, C2, C 3, C 4/4. 137 3, 1-042 3, -.20 84, -1.9 394/ 

DATA C5,C6 ,C7,C8,C9 /-2. 31 2 3, 1 . 7 187 *1 . 4 82 8, -2 . 785, 
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1-3.0368/ 

DATA C 10#C11»C12»C13/1.1318#2.586*«480*“. 19 6/ 

DATA C 14 *C l 5»C1 6>Ct 7/- 2.4243, U853444 5*-. 447 #-3.4 81/ 
DATA C02,C03*C04»C05,C06>C07/1.2930*0.2400#-.0982* 
1-.1762 ,.0607,.0494/C12P*C13P*C14P,C15P/-1.2121 

2».926?,.l651#".l 987/C22 *C23 *C24,C25/“1»7 4 06#-«2235 
Z* .2371#-. 1754/ 

DATA G AHA/ 1.2/ 

END 



